SyFi 0.3
|
00001 // Copyright (C) 2006-2009 Kent-Andre Mardal and Simula Research Laboratory. 00002 // Licensed under the GNU GPL Version 2, or (at your option) any later version. 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 // remove previously computed basis functions and dofs 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 // see e.g. Brezzi and Fortin book page 116 for the definition 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 // create the polynomial space 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 // create the dofs 00062 GiNaC::symbol t("t"); 00063 for (int j=0; j< 3; j++) 00064 { 00065 00066 // solve the linear system to compute 00067 // each of the basis functions 00068 GiNaC::lst equations; 00069 for (int i=0; i< 3; i++) 00070 { 00071 00072 Line line = triangle.line(i); 00073 // GiNaC::ex dofi = line.integrate(polynom); 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 // GiNaC::ex dofi = polynom.subs( x == midpoint.op(0)).subs( y == midpoint.op(1)); 00080 GiNaC::ex dofi = line.integrate(polynom); 00081 equations.append( dofi == dirac(i,j)); 00082 00083 if (j == 1) 00084 { 00085 // GiNaC::lst d = GiNaC::lst(line.vertex(0) , line.vertex(1)); 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 // dofs related to edges 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 // GiNaC::lst d = GiNaC::lst(triangle.vertex(0), triangle.vertex(1), triangle.vertex(2)); 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 // ------------VectorCrouzeixRaviart --- 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 } // namespace SyFi