SyFi 0.3
ArnoldFalkWintherWeakSym.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines