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 "ArnoldFalkWintherWeakSym.h" 00005 #include "Nedelec2Hdiv.h" 00006 #include "DiscontinuousLagrange.h" 00007 #include "P0.h" 00008 #include "utilities.h" 00009 00010 using std::cout; 00011 using std::endl; 00012 00013 namespace SyFi 00014 { 00015 00016 //------------ sigma element 00017 00018 ArnoldFalkWintherWeakSymSigma::ArnoldFalkWintherWeakSymSigma() : StandardFE() 00019 { 00020 description = "ArnoldFalkWintherWeakSymSigma"; 00021 } 00022 00023 ArnoldFalkWintherWeakSymSigma::ArnoldFalkWintherWeakSymSigma(Polygon& p, int order) : StandardFE(p, order) 00024 { 00025 compute_basis_functions(); 00026 } 00027 00028 void ArnoldFalkWintherWeakSymSigma:: compute_basis_functions() 00029 { 00030 00031 // remove previously computed basis functions and dofs 00032 Ns.clear(); 00033 dofs.clear(); 00034 00035 if ( order < 1 ) 00036 { 00037 throw(std::logic_error("Arnold-Falk-Winther elements must be of order 1 or higher.")); 00038 } 00039 00040 if ( p == NULL ) 00041 { 00042 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed")); 00043 } 00044 00045 description = istr("ArnoldFalkWintherWeakSymSigma_", order) + "_3D"; 00046 00047 Nedelec2Hdiv fe; 00048 fe.set_order(order); 00049 fe.set_polygon(*p); 00050 fe.compute_basis_functions(); 00051 00052 for (int d=0; d<3; d++) 00053 { 00054 for (unsigned int i=0; i<fe.nbf(); i++) 00055 { 00056 GiNaC::matrix Nmat = GiNaC::matrix(3,3); 00057 Nmat(d,0) = fe.N(i).op(0); 00058 Nmat(d,1) = fe.N(i).op(1); 00059 Nmat(d,2) = fe.N(i).op(2); 00060 Ns.insert(Ns.end(), Nmat); 00061 dofs.insert(dofs.end(),GiNaC::lst(fe.dof(i).op(0), fe.dof(i).op(1), d)); 00062 } 00063 } 00064 } 00065 00066 //------------ u element 00067 00068 ArnoldFalkWintherWeakSymU::ArnoldFalkWintherWeakSymU() : StandardFE() 00069 { 00070 description = "ArnoldFalkWintherWeakSymU"; 00071 } 00072 00073 ArnoldFalkWintherWeakSymU ::ArnoldFalkWintherWeakSymU (Polygon& p, int order) : StandardFE(p, order) 00074 { 00075 compute_basis_functions(); 00076 } 00077 00078 void ArnoldFalkWintherWeakSymU:: compute_basis_functions() 00079 { 00080 00081 // remove previously computed basis functions and dofs 00082 Ns.clear(); 00083 dofs.clear(); 00084 00085 if ( order < 1 ) 00086 { 00087 throw(std::logic_error("Arnold-Falk-Winther elements must be of order 1 or higher.")); 00088 } 00089 00090 if ( p == NULL ) 00091 { 00092 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed")); 00093 } 00094 00095 description = istr("ArnoldFalkWintherWeakSymU_", order) + "_3D"; 00096 00097 if ( order > 1 ) 00098 { 00099 VectorDiscontinuousLagrange fe; 00100 fe.set_order(order-1); 00101 fe.set_size(3); 00102 fe.set_polygon(*p); 00103 fe.compute_basis_functions(); 00104 00105 for (unsigned int i=0; i<fe.nbf(); i++) 00106 { 00107 GiNaC::lst Ni = GiNaC::lst(fe.N(i).op(0), fe.N(i).op(1), fe.N(i).op(2)); 00108 GiNaC::ex Nmat = GiNaC::matrix(3,1,Ni); 00109 Ns.insert(Ns.end(), Nmat); 00110 dofs.insert(dofs.end(),GiNaC::lst(fe.dof(i))); 00111 } 00112 } 00113 else if ( order == 1 ) 00114 { 00115 VectorP0 fe; 00116 fe.set_order(order-1); 00117 fe.set_size(3); 00118 fe.set_polygon(*p); 00119 fe.compute_basis_functions(); 00120 00121 for (unsigned int i=0; i<fe.nbf(); i++) 00122 { 00123 GiNaC::lst Ni = GiNaC::lst(fe.N(i).op(0), fe.N(i).op(1), fe.N(i).op(2)); 00124 GiNaC::ex Nmat = GiNaC::matrix(3,1,Ni); 00125 Ns.insert(Ns.end(), Nmat); 00126 GiNaC::ex d = GiNaC::lst(fe.dof(i), 0); 00127 dofs.insert(dofs.end(),GiNaC::lst(d)); 00128 } 00129 } 00130 00131 } 00132 00133 //------------ p element 00134 00135 ArnoldFalkWintherWeakSymP::ArnoldFalkWintherWeakSymP() : StandardFE() 00136 { 00137 description = "ArnoldFalkWintherWeakSymP"; 00138 } 00139 00140 ArnoldFalkWintherWeakSymP ::ArnoldFalkWintherWeakSymP (Polygon& p, int order) : StandardFE(p, order) 00141 { 00142 compute_basis_functions(); 00143 } 00144 00145 void ArnoldFalkWintherWeakSymP:: compute_basis_functions() 00146 { 00147 00148 // remove previously computed basis functions and dofs 00149 Ns.clear(); 00150 dofs.clear(); 00151 00152 if ( order < 1 ) 00153 { 00154 cout <<"Arnold-Falk-Winther elements must be of order 1 or higher."<<endl; 00155 return; 00156 } 00157 00158 if ( p == NULL ) 00159 { 00160 cout <<"You need to set a polygon before the basisfunctions can be computed"<<endl; 00161 return; 00162 } 00163 00164 description = istr("ArnoldFalkWintherWeakSymP_", order) + "_3D"; 00165 00166 if ( order > 1 ) 00167 { 00168 00169 VectorDiscontinuousLagrange fe; 00170 fe.set_order(order); 00171 fe.set_size(3); 00172 fe.set_polygon(*p); 00173 fe.compute_basis_functions(); 00174 00175 for (unsigned int i=0; i<fe.nbf(); i++) 00176 { 00177 GiNaC::matrix Nmat = GiNaC::matrix(3,3); 00178 Nmat(1,2) = -fe.N(i).op(0); 00179 Nmat(2,1) = fe.N(i).op(0); 00180 00181 Nmat(0,2) = fe.N(i).op(1); 00182 Nmat(2,0) = -fe.N(i).op(1); 00183 00184 Nmat(0,1) = -fe.N(i).op(2); 00185 Nmat(1,0) = fe.N(i).op(2); 00186 00187 Ns.insert(Ns.end(), Nmat); 00188 dofs.insert(dofs.end(),GiNaC::lst(fe.dof(i))); 00189 } 00190 } 00191 else if ( order == 1 ) 00192 { 00193 00194 VectorP0 fe; 00195 fe.set_order(order); 00196 fe.set_size(3); 00197 fe.set_polygon(*p); 00198 fe.compute_basis_functions(); 00199 00200 for (unsigned int i=0; i<fe.nbf(); i++) 00201 { 00202 GiNaC::matrix Nmat = GiNaC::matrix(3,3); 00203 Nmat(1,2) = -fe.N(i).op(0); 00204 Nmat(2,1) = fe.N(i).op(0); 00205 00206 Nmat(0,2) = fe.N(i).op(1); 00207 Nmat(2,0) = -fe.N(i).op(1); 00208 00209 Nmat(0,1) = -fe.N(i).op(2); 00210 Nmat(1,0) = fe.N(i).op(2); 00211 00212 Ns.insert(Ns.end(), Nmat); 00213 GiNaC::ex d = GiNaC::lst(fe.dof(i), 0); 00214 dofs.insert(dofs.end(),GiNaC::lst(d)); 00215 } 00216 } 00217 00218 } 00219 00220 } // namespace SyFi