arnoldfalkwinther.cpp

Go to the documentation of this file.

#include <SyFi.h>

using namespace GiNaC; 

int main() {
   ReferenceTetrahedron tetrahedron; 

   nsd = 3; 
   int order = 3; 


   ArnoldFalkWintherWeakSymSigma sigma_fe; 
   sigma_fe.set(order); 
   sigma_fe.set(tetrahedron); 
   sigma_fe.compute_basis_functions(); 
   for (int i=0; i<sigma_fe.nbf(); i++) {
     cout <<"sigma_fe.N("<<i<<")="<<sigma_fe.N(i)<<endl; ; 
     cout <<"div(sigma_fe.N("<<i<<"))="<<div(sigma_fe.N(i))<<endl; ; 
     cout <<"sigma_fe.dof("<<i<<"))="<<sigma_fe.dof(i)<<endl; ; 
   }

   map<pair<int,int>,ex> A; 
   pair<int,int> index; 
   for (int i=0; i<sigma_fe.nbf(); i++) {
     index.first = i; 
     for (int j=0; j<sigma_fe.nbf(); j++) {
       index.second = j; 
       ex integrand = inner(sigma_fe.N(i), sigma_fe.N(j)); 
       A[index] = tetrahedron.integrate(integrand); 
       cout <<"A["<<i<<","<<j<<"]="<<A[index]<<endl; 
     }
   }




   ArnoldFalkWintherWeakSymU u_fe; 
   u_fe.set(order); 
   u_fe.set(tetrahedron); 
   u_fe.compute_basis_functions(); 
   for (int i=0; i<u_fe.nbf(); i++) {
     cout <<"u_fe.N("<<i<<")="<<u_fe.N(i)<<endl; ; 
   }

   map<pair<int,int>,ex> B; 
   for (int i=0; i<sigma_fe.nbf(); i++) {
     index.first = i; 
     for (int j=0; j<u_fe.nbf(); j++) {
       index.second = j; 
       ex integrand = inner(div(sigma_fe.N(i)), u_fe.N(j)); 
       B[index] = tetrahedron.integrate(integrand); 
       cout <<"B["<<i<<","<<j<<"]="<<B[index]<<endl; 
     }
   }




   ArnoldFalkWintherWeakSymP p_fe; 
   p_fe.set(order); 
   p_fe.set(tetrahedron); 
   p_fe.compute_basis_functions(); 
   for (int i=0; i<p_fe.nbf(); i++) {
     cout <<"p_fe.N("<<i<<")="<<p_fe.N(i)<<endl; ; 
   }



}

---