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 "P0.h" 00005 #include "tools.h" 00006 00007 using std::cout; 00008 using std::endl; 00009 00010 namespace SyFi 00011 { 00012 00013 P0::P0() : StandardFE() 00014 { 00015 description = "P0"; 00016 } 00017 00018 P0:: P0(Polygon& p, unsigned int order) : StandardFE (p,order) 00019 { 00020 compute_basis_functions(); 00021 } 00022 00023 void P0:: compute_basis_functions() 00024 { 00025 00026 // remove previously computed basis functions and dofs 00027 Ns.clear(); 00028 dofs.clear(); 00029 00030 if ( p == NULL ) 00031 { 00032 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed")); 00033 } 00034 00035 // insert basis function 00036 Ns.insert(Ns.end(), GiNaC::numeric(1)); 00037 00038 GiNaC::lst midpoint = GiNaC::lst(); 00039 // create and insert dof 00040 // p is a lst 00041 if (GiNaC::is_a<GiNaC::lst>(p->vertex(0))) 00042 { 00043 for (unsigned int d=0; d< p->vertex(1).nops(); d++) 00044 { 00045 midpoint.append(GiNaC::numeric(0)); 00046 } 00047 for (unsigned int i=0; i< p->no_vertices(); i++) 00048 { 00049 int nops; 00050 nops = p->vertex(i).nops(); 00051 for (int d=0; d< nops; d++) 00052 { 00053 midpoint.let_op(d) += p->vertex(i).op(d); 00054 } 00055 } 00056 for (unsigned int d=0; d< p->vertex(1).nops(); d++) 00057 { 00058 midpoint.let_op(d) = midpoint.op(d)/p->no_vertices(); 00059 } 00060 } 00061 else 00062 { 00063 midpoint.append(GiNaC::numeric(0)); 00064 for (unsigned int i=0; i< p->no_vertices(); i++) 00065 { 00066 midpoint.let_op(0) += p->vertex(i); 00067 } 00068 midpoint.let_op(0) = midpoint.op(0)/p->no_vertices(); 00069 } 00070 00071 dofs.insert(dofs.end(), midpoint); 00072 00073 description = istr("P0_" , midpoint.nops()) + "D"; 00074 } 00075 00076 // ------------VectorP0 --- 00077 00078 VectorP0::VectorP0() : StandardFE() 00079 { 00080 description = "VectorP0"; 00081 } 00082 00083 VectorP0::VectorP0(Polygon& p, unsigned int order, unsigned int size_) : StandardFE(p, order) 00084 { 00085 size = size_ < 0 ? nsd: size_; 00086 compute_basis_functions(); 00087 } 00088 00089 void VectorP0:: compute_basis_functions() 00090 { 00091 00092 // remove previously computed basis functions and dofs 00093 Ns.clear(); 00094 dofs.clear(); 00095 00096 if ( p == NULL ) 00097 { 00098 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed")); 00099 } 00100 00101 if ( size == 0) 00102 { 00103 throw(std::logic_error("You need to set the size of the vector before the basisfunctions can be computed")); 00104 } 00105 00106 P0 fe(*p); 00107 GiNaC::lst zero_list; 00108 for (unsigned int s=1; s<= size ; s++) 00109 { 00110 zero_list.append(0); 00111 } 00112 00113 for (unsigned int i=0; i< fe.nbf() ; i++) 00114 { 00115 for (unsigned int s=0; s< size ; s++) 00116 { 00117 GiNaC::lst Nis = zero_list; 00118 Nis.let_op(s) = fe.N(i); 00119 GiNaC::ex Nmat = GiNaC::matrix(size,1,Nis); 00120 Ns.insert(Ns.end(), Nmat); 00121 00122 GiNaC::lst dof = GiNaC::lst(fe.dof(i), s) ; 00123 dofs.insert(dofs.end(), dof); 00124 } 00125 } 00126 description = "Vector" + fe.str(); 00127 } 00128 00129 void VectorP0:: set_size(unsigned int size_) 00130 { 00131 size = size_; 00132 } 00133 00134 // ------------TensorP0 --- 00135 00136 TensorP0::TensorP0() : StandardFE() 00137 { 00138 description = "TensorP0"; 00139 } 00140 00141 TensorP0::TensorP0(Polygon& p, unsigned int order, unsigned int size_) : StandardFE(p, order) 00142 { 00143 size = size_ < 0 ? nsd: size_; 00144 compute_basis_functions(); 00145 } 00146 00147 void TensorP0:: compute_basis_functions() 00148 { 00149 00150 // remove previously computed basis functions and dofs 00151 Ns.clear(); 00152 dofs.clear(); 00153 00154 if ( p == NULL ) 00155 { 00156 throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed")); 00157 } 00158 00159 if ( size == 0) 00160 { 00161 throw(std::logic_error("You need to set the size of the vector before the basisfunctions can be computed")); 00162 } 00163 00164 P0 fe(*p); 00165 GiNaC::lst zero_list; 00166 for (unsigned int s=1; s<= size*size ; s++) 00167 { 00168 zero_list.append(0); 00169 } 00170 00171 for (unsigned int i=0; i< fe.nbf() ; i++) 00172 { 00173 for (unsigned int r=0; r< size ; r++) 00174 { 00175 for (unsigned int s=0; s< size ; s++) 00176 { 00177 GiNaC::lst Nis = zero_list; 00178 Nis.let_op((size)*r + s) = fe.N(i); 00179 GiNaC::ex Nmat = GiNaC::matrix(size,size,Nis); 00180 Ns.insert(Ns.end(), Nmat); 00181 00182 GiNaC::lst dof = GiNaC::lst(fe.dof(i), r, s) ; 00183 dofs.insert(dofs.end(), dof); 00184 } 00185 } 00186 } 00187 description = "Tensor" + fe.str(); 00188 } 00189 00190 void TensorP0:: set_size(unsigned int size_) 00191 { 00192 size = size_; 00193 } 00194 00195 }