CrouzeixRaviart.cpp
Go to the documentation of this file.00001
00002
00003
00004 #include "CrouzeixRaviart.h"
00005
00006 #include "tools.h"
00007
00008 using std::cout;
00009 using std::endl;
00010 using std::string;
00011
00012 namespace SyFi
00013 {
00014
00015 CrouzeixRaviart:: CrouzeixRaviart() : StandardFE()
00016 {
00017 order = 1;
00018 }
00019
00020 CrouzeixRaviart:: CrouzeixRaviart (Polygon& p, unsigned int order) : StandardFE(p, order)
00021 {
00022 compute_basis_functions();
00023 }
00024
00025 void CrouzeixRaviart:: compute_basis_functions()
00026 {
00027
00028
00029 Ns.clear();
00030 dofs.clear();
00031
00032 if ( p == NULL )
00033 {
00034 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
00035 }
00036
00037 if (order != 1)
00038 {
00039 throw(std::logic_error("Only Crouziex-Raviart elements of order 1 is possible"));
00040 }
00041
00042
00043
00044 if ( p->str().find("ReferenceLine") != string::npos )
00045 {
00046 cout <<"Can not define the Raviart-Thomas element on a line"<<endl;
00047 }
00048 else if ( p->str().find("Triangle") != string::npos )
00049 {
00050
00051 description = "CrouzeixRaviart_2D";
00052
00053 Triangle& triangle = (Triangle&)(*p);
00054
00055
00056 GiNaC::ex polynom_space = bernstein(1, triangle, "a");
00057 GiNaC::ex polynom = polynom_space.op(0);
00058 GiNaC::lst variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00059 GiNaC::ex basis = polynom_space.op(2);
00060
00061
00062 GiNaC::symbol t("t");
00063 for (int j=0; j< 3; j++)
00064 {
00065
00066
00067
00068 GiNaC::lst equations;
00069 for (int i=0; i< 3; i++)
00070 {
00071
00072 Line line = triangle.line(i);
00073
00074 GiNaC::lst midpoint = GiNaC::lst(
00075 (line.vertex(0).op(0) + line.vertex(1).op(0))/2,
00076 (line.vertex(0).op(1) + line.vertex(1).op(1))/2);
00077 dofs.insert(dofs.end(), midpoint);
00078
00079
00080 GiNaC::ex dofi = line.integrate(polynom);
00081 equations.append( dofi == dirac(i,j));
00082
00083 if (j == 1)
00084 {
00085
00086 dofs.insert(dofs.end(), midpoint);
00087 }
00088
00089 }
00090 GiNaC::ex sub = lsolve(equations, variables);
00091 GiNaC::ex Ni = polynom.subs(sub);
00092 Ns.insert(Ns.end(),Ni);
00093 }
00094
00095 }
00096 else if ( p->str().find("Tetrahedron") != string::npos )
00097 {
00098
00099 description = "CrouzeixRaviart_3D";
00100
00101 Tetrahedron& tetrahedron = (Tetrahedron&)(*p);
00102 GiNaC::ex polynom_space = bernstein(1, tetrahedron, "a");
00103 GiNaC::ex polynom = polynom_space.op(0);
00104 GiNaC::lst variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
00105 GiNaC::ex basis = polynom_space.op(2);
00106
00107 GiNaC::ex bernstein_pol;
00108
00109 GiNaC::symbol t("t");
00110
00111 for (int j=0; j< 4; j++)
00112 {
00113
00114 GiNaC::lst equations;
00115 for (int i=0; i< 4; i++)
00116 {
00117 Triangle triangle = tetrahedron.triangle(i);
00118 GiNaC::lst midpoint = GiNaC::lst(
00119 (triangle.vertex(0).op(0) + triangle.vertex(1).op(0) + triangle.vertex(2).op(0))/3,
00120 (triangle.vertex(0).op(1) + triangle.vertex(1).op(1) + triangle.vertex(2).op(1))/3,
00121 (triangle.vertex(0).op(2) + triangle.vertex(1).op(2) + triangle.vertex(2).op(2))/3
00122 );
00123
00124 GiNaC::ex dofi = polynom.subs(x == midpoint.op(0)).subs(y == midpoint.op(1)).subs(z == midpoint.op(2));
00125 equations.append( dofi == dirac(i,j));
00126
00127 if ( j == 1 )
00128 {
00129
00130 dofs.insert(dofs.end(), midpoint);
00131 }
00132 }
00133 GiNaC::ex sub = lsolve(equations, variables);
00134 GiNaC::ex Ni = polynom.subs(sub);
00135 Ns.insert(Ns.end(),Ni);
00136
00137 }
00138 }
00139 }
00140
00141
00142
00143 VectorCrouzeixRaviart:: VectorCrouzeixRaviart (Polygon& p, unsigned int order, unsigned int size_) : StandardFE(p, order)
00144 {
00145 size = size_ < 0 ? nsd: size_;
00146 compute_basis_functions();
00147 }
00148
00149 VectorCrouzeixRaviart:: VectorCrouzeixRaviart() : StandardFE()
00150 {
00151 description = "VectorCrouzeixRaviart";
00152 order = 1;
00153 }
00154
00155 void VectorCrouzeixRaviart:: compute_basis_functions()
00156 {
00157
00158 if (order != 1)
00159 {
00160 throw(std::logic_error("Only Crouziex-Raviart elements of order 1 is possible"));
00161 }
00162
00163 CrouzeixRaviart fe;
00164 fe.set_polygon(*p);
00165 fe.compute_basis_functions();
00166
00167 description = "Vector" + fe.str();
00168
00169 GiNaC::lst zero_list;
00170 for (unsigned int s=1; s<= size ; s++)
00171 {
00172 zero_list.append(0);
00173 }
00174
00175 for (unsigned int s=0; s< size ; s++)
00176 {
00177 for (unsigned int i=0; i< fe.nbf() ; i++)
00178 {
00179 GiNaC::lst Nis = zero_list;
00180 Nis.let_op(s) = fe.N(i);
00181 GiNaC::ex Nmat = GiNaC::matrix(size,1,Nis);
00182 Ns.insert(Ns.end(), Nmat);
00183
00184 GiNaC::lst dof = GiNaC::lst(fe.dof(i), s) ;
00185 dofs.insert(dofs.end(), dof);
00186 }
00187 }
00188 }
00189
00190 void VectorCrouzeixRaviart:: set_size(unsigned int size_)
00191 {
00192 size = size_;
00193 }
00194
00195 }