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 "ElementComputations.h" 00005 #include "tools.h" 00006 00007 using std::cout; 00008 using std::endl; 00009 00010 namespace SyFi 00011 { 00012 00013 void usage(FE& fe) 00014 { 00015 for (unsigned int i=0; i< fe.nbf(); i++) 00016 { 00017 cout <<"fe.N("<<i<<") = "<<fe.N(i)<<endl; 00018 cout <<"grad(fe.N("<<i<<")) = "<<grad(fe.N(i))<<endl; 00019 cout <<"fe.dof("<<i<<") = "<<fe.dof(i)<<endl; 00020 } 00021 } 00022 00023 void usage(FE& v_fe, FE& p_fe) 00024 { 00025 for (unsigned int i=0; i< v_fe.nbf(); i++) 00026 { 00027 cout <<"v_fe.N("<<i<<") = "<<v_fe.N(i)<<endl; 00028 cout <<"grad(v_fe.N("<<i<<")) = "<<grad(v_fe.N(i))<<endl; 00029 cout <<"v_fe.dof("<<i<<") = "<<v_fe.dof(i)<<endl; 00030 } 00031 for (unsigned int i=0; i< p_fe.nbf(); i++) 00032 { 00033 cout <<"p_fe.N("<<i<<")= "<<p_fe.N(i)<<endl; 00034 cout <<"p_fe.dof("<<i<<")= "<<p_fe.dof(i)<<endl; 00035 } 00036 } 00037 00038 void compute_Poisson_element_matrix( 00039 FE& fe, 00040 Dof& dof, 00041 std::map<std::pair<unsigned int,unsigned int>, GiNaC::ex>& A) 00042 { 00043 std::pair<unsigned int,unsigned int> index; 00044 00045 // Insert the local degrees of freedom into the global Dof 00046 for (unsigned int i=0; i< fe.nbf(); i++) 00047 { 00048 dof.insert_dof(1,i,fe.dof(i)); 00049 } 00050 00051 Polygon& domain = fe.get_polygon(); 00052 00053 // The term (grad u, grad v) 00054 for (unsigned int i=0; i< fe.nbf(); i++) 00055 { 00056 // fetch the global dof for Ni 00057 index.first = dof.glob_dof(fe.dof(i)); 00058 for (unsigned int j=0; j< fe.nbf(); j++) 00059 { 00060 // fetch the global dof for Nj 00061 index.second = dof.glob_dof(fe.dof(j)); 00062 // compute the integrand 00063 GiNaC::ex nabla = inner(grad(fe.N(i)), 00064 grad(fe.N(j))); 00065 // compute the integral 00066 GiNaC::ex Aij = domain.integrate(nabla); 00067 A[index] += Aij; // add to global matrix 00068 } 00069 } 00070 } 00071 00072 void compute_Stokes_element_matrix( 00073 FE& v_fe, 00074 FE& p_fe, 00075 Dof& dof, 00076 std::map<std::pair<unsigned int,unsigned int>, GiNaC::ex>& A) 00077 { 00078 std::pair<unsigned int,unsigned int> index; 00079 std::pair<unsigned int,unsigned int> index2; 00080 00081 // FIXME: need to check that p_fe 00082 // contains the same domain 00083 Polygon& domain = v_fe.get_polygon(); 00084 00085 // Insert the local degrees of freedom into the global Dof 00086 for (unsigned int i=0; i< v_fe.nbf(); i++) 00087 { 00088 dof.insert_dof(1,i,v_fe.dof(i)); 00089 } 00090 for (unsigned int i=0; i< p_fe.nbf(); i++) 00091 { 00092 dof.insert_dof(1,v_fe.nbf()+i,p_fe.dof(i)); 00093 } 00094 00095 // The term (grad u, grad v) 00096 for (unsigned int i=0; i< v_fe.nbf(); i++) 00097 { 00098 // fetch the global dof for v_i 00099 index.first = dof.glob_dof(v_fe.dof(i)); 00100 for (unsigned int j=0; j< v_fe.nbf(); j++) 00101 { 00102 // fetch the global dof for v_j 00103 index.second = dof.glob_dof(v_fe.dof(j)); 00104 GiNaC::ex nabla = inner(grad(v_fe.N(i)), 00105 grad(v_fe.N(j)));// compute the integrand 00106 // compute the integral 00107 GiNaC::ex Aij = domain.integrate(nabla); 00108 A[index] += Aij; // add to global matrix 00109 } 00110 } 00111 00112 // The term -(div u, q) 00113 for (unsigned int i=0; i< p_fe.nbf(); i++) 00114 { 00115 // fetch the global dof for p_i 00116 index.first = dof.glob_dof(p_fe.dof(i)); 00117 for (unsigned int j=0; j< v_fe.nbf(); j++) 00118 { 00119 // fetch the global dof for v_j 00120 index.second=dof.glob_dof(v_fe.dof(j)); 00121 // compute the integrand 00122 GiNaC::ex divV= -p_fe.N(i)*div(v_fe.N(j)); 00123 // compute the integral 00124 GiNaC::ex Aij = domain.integrate(divV); 00125 A[index] += Aij; // add to global matrix 00126 00127 // Do not need to compute the term (grad(p),v), since the system is 00128 // symmetric. We simply set Aji = Aij 00129 index2.first = index.second; 00130 index2.second = index.first; 00131 A[index2] += Aij; 00132 } 00133 } 00134 } 00135 00136 void compute_mixed_Poisson_element_matrix( 00137 FE& v_fe, 00138 FE& p_fe, 00139 Dof& dof, 00140 std::map<std::pair<unsigned int,unsigned int>, GiNaC::ex>& A) 00141 { 00142 std::pair<unsigned int,unsigned int> index; 00143 std::pair<unsigned int,unsigned int> index2; 00144 00145 // FIXME: need to check that p_fe 00146 // contains the same domain 00147 Polygon& domain = v_fe.get_polygon(); 00148 00149 // Insert the local degrees of freedom into the global Dof 00150 for (unsigned int i=0; i< v_fe.nbf(); i++) 00151 { 00152 dof.insert_dof(1,i,v_fe.dof(i)); 00153 } 00154 for (unsigned int i=0; i< p_fe.nbf(); i++) 00155 { 00156 dof.insert_dof(1,v_fe.nbf()+i+1,p_fe.dof(i)); 00157 } 00158 00159 // The term (u,v) 00160 for (unsigned int i=0; i< v_fe.nbf(); i++) 00161 { 00162 // fetch the global dof related to i and v 00163 index.first = dof.glob_dof(v_fe.dof(i)); 00164 for (unsigned int j=0; j< v_fe.nbf(); j++) 00165 { 00166 // fetch the global dof related to j and p 00167 index.second = dof.glob_dof(v_fe.dof(j)); 00168 // compute the integrand 00169 GiNaC::ex mass = inner(v_fe.N(i),v_fe.N(j)); 00170 // compute the integral 00171 GiNaC::ex Aij = domain.integrate(mass); 00172 A[index] += Aij; // add to global matrix 00173 } 00174 } 00175 00176 // The term -(div u, q) 00177 for (unsigned int i=0; i< p_fe.nbf(); i++) 00178 { 00179 // fetch the global dof for p_i 00180 index.first = dof.glob_dof(p_fe.dof(i)); 00181 for (unsigned int j=0; j< v_fe.nbf(); j++) 00182 { 00183 // fetch the global dof for v_j 00184 index.second=dof.glob_dof(v_fe.dof(j)); 00185 // compute the integrand 00186 GiNaC::ex divV= -p_fe.N(i)*div(v_fe.N(j)); 00187 // compute the integral 00188 GiNaC::ex Aij = domain.integrate(divV); 00189 A[index] += Aij; // add to global matrix 00190 00191 // Do not need to compute the term (grad(p),v), since the system is 00192 // symmetric we simply set Aji = Aij 00193 index2.first = index.second; 00194 index2.second = index.first; 00195 A[index2] += Aij; 00196 } 00197 } 00198 } 00199 00200 }