-
Notifications
You must be signed in to change notification settings - Fork 7
/
Poisson.cpp
153 lines (136 loc) · 5.28 KB
/
Poisson.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
#include <iostream>
#include "ProjectFEA.hpp"
/// A new Model is defined here. This weak form is integrated over each element.
/// The Virtual 'weak_form' function defined in 'LocalIntegrator' is overloaded during runtime.
class new_LocalIntegrator: public LocalIntegrator<TrialFunction>
{
public:
new_LocalIntegrator(Form<TrialFunction>& a, TrialFunction& u, TestFunctionGalerkin<TrialFunction>& v):
LocalIntegrator (a,u,v)
{
}
sp_mat weak_form(Form<TrialFunction>& a, TrialFunction& u, TestFunctionGalerkin<TrialFunction>& v)
{
return a.inner(a.grad(v),a.grad(u))*a.dX(u);
}
};
/// The Virtual 'weak_form_vector' function defined in 'LocalIntegrator' is overloaded during runtime.
class new_LocalIntegrator2: public LocalIntegrator<TrialFunction>
{
public:
new_LocalIntegrator2(Form<TrialFunction>& a, TrialFunction& u, TestFunctionGalerkin<TrialFunction>& v):
LocalIntegrator (a,u,v)
{
}
mat weak_form_vector(Form<TrialFunction>& a, TrialFunction& u, TestFunctionGalerkin<TrialFunction>& v)
{
vec b;
b<<0<<endr<<0<<endr<<0<<endr;
return a.dot(v,b)*a.dX(u);
}
};
/// The Virtual 'weak_form_vector' function defined in 'LocalIntegrator' is overloaded during runtime.
class new_Neu_Surf_LclIntgrtr: public LocalIntegrator <TrialFunctionNeumannSurface>
{
public:
new_Neu_Surf_LclIntgrtr(Form<TrialFunctionNeumannSurface>& a, TrialFunctionNeumannSurface& u,
TestFunctionGalerkin<TrialFunctionNeumannSurface>& v):
LocalIntegrator (a,u,v)
{
}
mat weak_form_vector(Form<TrialFunctionNeumannSurface>& a, TrialFunctionNeumannSurface& u,
TestFunctionGalerkin<TrialFunctionNeumannSurface>& v)
{
vec vctr={1,1,1};
return a.dot(v,vctr)*a.dS(u);
}
};
/// The Virtual 'weak_form_vector' function defined in 'LocalIntegrator' is overloaded during runtime.
class new_Neu_Line_LclIntgrtr: public LocalIntegrator<TrialFunctionNeumannLine>
{
public:
new_Neu_Line_LclIntgrtr(Form<TrialFunctionNeumannLine>& a, TrialFunctionNeumannLine& u,
TestFunctionGalerkin<TrialFunctionNeumannLine>& v):
LocalIntegrator (a,u,v)
{
}
mat weak_form_vector(Form<TrialFunctionNeumannLine>& a, TrialFunctionNeumannLine& u,
TestFunctionGalerkin<TrialFunctionNeumannLine>& v)
{
//vec vctr=vec(a.x(u));
vec vctr={1,1,1};
return a.dot(v,vctr)*a.dL(u);
}
};
///This class over here through its overloaded virtual function declares the values of Dirichlet Nodes.
/// The virtual function 'Eval' is evaluated at each node to find the value of Dirichlet Condtion at that node.
class DeclaredExprssn : public Expression
{
public:
DeclaredExprssn (int vectorLevel): Expression (vectorLevel)
{
}
vec Eval(vec& x)
{
vec value={0,0,0};
return value;
}
};
int main(int argc, char *argv[])
{
if(argc<3 || argc>3)
{
std::cout<<"Usage: ./Poisson <.msh Filename> <Dimension>\n";
return 0;
}
std::string FileName(argv[1]);
int Dimension=*argv[2]-'0';
//cout<<"File Name= "<<FileName<<"\n";
//cout<<"Dimension= "<<Dimension<<"\n";
libGmshReader::MeshReader Mesh(FileName, Dimension);
int vectorLevel=3;
/// Declaration of Quantity to be Calculated.
TrialFunction u(Mesh,vectorLevel);
TestFunctionGalerkin<TrialFunction> v(u);
Form<TrialFunction> a;
new_LocalIntegrator lcl_intgrt(a,u,v);
SystemAssembler<new_LocalIntegrator, TrialFunction> systmAssmbly(a,u,v);
//systmAssmbly.SetLocalIntegrator(lcl_intgrt);
//sp_mat A;
VariableMatrix A(1);
systmAssmbly.SetMatrixSize(A.Matrix[0][0]);
systmAssmbly.RunSystemAssembly(lcl_intgrt, A.Matrix[0][0]);
Form<TrialFunction> a2;
new_LocalIntegrator2 lcl_intgrt2(a2,u,v);
SystemAssembler<new_LocalIntegrator2, TrialFunction> systmAssmbly2(a2,u,v);
//mat b;
VariableVector b(1);
systmAssmbly2.SetVectorSize(b.Vector[0]);
systmAssmbly2.RunSystemAssemblyVector(lcl_intgrt2,b.Vector[0]);
Form<TrialFunctionNeumannSurface> a3;
TrialFunctionNeumannSurface u_surf(u,0);
TestFunctionGalerkin<TrialFunctionNeumannSurface> v_surf(u_surf);
new_Neu_Surf_LclIntgrtr lclintgtr3(a3,u_surf, v_surf);
SystemAssembler<new_Neu_Surf_LclIntgrtr, TrialFunctionNeumannSurface> systmAssmbly3(a3, u_surf, v_surf);
systmAssmbly3.RunSystemAssemblyVector(lclintgtr3, b.Vector[0]);
//Form<TrialFunctionNeumannLine> a4;
//TrialFunctionNeumannLine u_line(u,0);
/*TestFunctionGalerkin<TrialFunctionNeumannLine> v_line(u_line);
new_Neu_Line_LclIntgrtr lcl_intgrt4(a4,u_line,v_line);
SystemAssembler<new_Neu_Line_LclIntgrtr, TrialFunctionNeumannLine> systmAssmbly4(a4,u_line, v_line);
systmAssmbly4.RunSystemAssemblyVector(lcl_intgrt4,b.Vector[0]);*/
//cout<<b.Vector[0];
umat boolDiricletNodes={1,1,1};
DirichletBC DrchltBC(u_surf,1, boolDiricletNodes);
DeclaredExprssn Dirich(vectorLevel);
DrchltBC.SetDirichletBCExpression(Dirich);
DrchltBC.ApplyBC(A.Matrix[0][0],b.Vector[0]);
//cout<<b.Vector[0];
mat X=spsolve(A.Matrix[0][0],b.Vector[0]);
//cout<<X;
GmshWriter Write(u, "poisson.pos");
Write.WriteToGmsh(X);
//Write.WriteToGmsh(X,"Disp");
cout<<"Done!!!\n";
return 0;
}