00001
00002
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
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
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
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
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
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
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 }