ElementComputations.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 "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 }

Generated on Mon Aug 31 16:16:45 2009 for SyFi by  doxygen 1.5.9