/[ff3d]/ff3d/solver/FEMDiscretization.hpp
ViewVC logotype

Diff of /ff3d/solver/FEMDiscretization.hpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.4 by delpinux, Sat Apr 12 17:09:16 2003 UTC revision 1.5 by delpinux, Tue Apr 15 17:47:27 2003 UTC
# Line 17  Line 17 
17    
18  //  $Id$  //  $Id$
19    
   
20  #ifndef FEM_DISCRETIZATION_HPP  #ifndef FEM_DISCRETIZATION_HPP
21  #define FEM_DISCRETIZATION_HPP  #define FEM_DISCRETIZATION_HPP
22    
# Line 40  Line 39 
39    
40  #include <ElementaryMatrixSet.hpp>  #include <ElementaryMatrixSet.hpp>
41    
42  #include <Discretization.hpp>  #include <BaseFEMDiscretization.hpp>
43    
44  #include <Mesh.hpp>  #include <Mesh.hpp>
45  #include <Structured3DMesh.hpp>  #include <Structured3DMesh.hpp>
46    
 #include <DiscretizedOperators.hpp>  
 #include <DegreeOfFreedomSet.hpp>  
   
 #include <FiniteElementTraits.hpp>  
   
 #include <PDE.hpp>  
 #include <PDEProblem.hpp>  
 #include <MassOperator.hpp>  
 #include <FirstOrderOperator.hpp>  
 #include <DivMuGrad.hpp>  
 #include <SecondOrderOperator.hpp>  
   
47  #include <Timer.hpp>  #include <Timer.hpp>
48  #include <Q1FiniteElement.hpp>  #include <Q1FiniteElement.hpp>
49    
50  #include <DoubleHashedMatrix.hpp>  #include <DoubleHashedMatrix.hpp>
51  #include <UnAssembledMatrix.hpp>  #include <UnAssembledMatrix.hpp>
52    
 #include <VariationalProblem.hpp>  
   
 #include <VariationalOperatorFV.hpp>  
 #include <VariationalOperatorFdxGV.hpp>  
 #include <VariationalOperatorFgradGgradV.hpp>  
   
 #include <ConformTransformation.hpp>  
   
53  #warning Should not use language classes here  #warning Should not use language classes here
54  #include <FunctionExpression.hpp>  #include <FunctionExpression.hpp>
55  #include <MeshExpression.hpp>  #include <MeshExpression.hpp>
# Line 78  Line 57 
57    
58  template <typename GivenMeshType>  template <typename GivenMeshType>
59  class FEMDiscretization  class FEMDiscretization
60    : public Discretization    : public BaseFEMDiscretization<GivenMeshType>
61  {  {
62  public:  private:
63    /// The type of mesh used for discretization    /// The type of mesh used for discretization
64    typedef GivenMeshType MeshType;    typedef GivenMeshType MeshType;
65    
# Line 105  public: Line 84  public:
84    /// Associated jacobian    /// Associated jacobian
85    typedef typename FiniteElement::JacobianTransformation   JacobianTransformation;    typedef typename FiniteElement::JacobianTransformation   JacobianTransformation;
86    
 private:  
   /// Mesh used to perform discretization  
   MeshType& __mesh;  
   
   /// Set of elementary matrices  
   mutable ElementaryMatrixSet <ElementaryMatrixType> __eSet;  
   
   /// Operators that are discretized  
   mutable DiscretizedOperators<ElementaryMatrixType> __discretizedOperators;  
   
   /// Set of degrees of freedom  
   const DegreeOfFreedomSet& __degreeOfFreedomSet;  
   
   /**  
    * Generates elementary vector  
    *  
    * @param eVector the generated elementary vector  
    * @param J the jacobian of the transformation  
    * @param f the function to discretize  
    */  
   void  
   generatesElementaryVector(ElementaryVectorType& eVector,  
                             const JacobianTransformation& J,  
                             const ElementaryVectorType& f) const  
   {  
     FiniteElementType::instance().integrateWj(eVector,J,f);  
   }  
   
   /**  
    * Generates elementary matrices set for a given element  
    *  
    * @param eSet the set of elementary matrices  
    * @param J the jacobian of the transformation  
    */  
   void  
   generatesElementaryMatrix(ElementaryMatrixSet<ElementaryMatrixType>& eSet,  
                             const JacobianTransformation& J) const  
   {  
     if (eSet.isMassOperator()) {  
       generatesElementaryMatrix(PDEOperator::massop,  
                                 J, eSet.massOperator());  
     }  
     if (eSet.isFirstOrderOperator()) {  
       for (size_t i=0; i<3; ++i) {  
         if (eSet.isFirstOrderUdxV(i)) {  
           generatesElementaryMatrix(PDEOperator::firstorderopTransposed,  
                                     J,eSet.firstOrderOperatorUdxV(i),i);  
         }  
         if (eSet.isFirstOrderDxUV(i)) {  
           generatesElementaryMatrix(PDEOperator::firstorderop, J,  
                                     eSet.firstOrderOperatorDxUV(i),i);  
         }  
       }  
     }  
     if (eSet.isSecondOrderOperator()) {  
       for (size_t i=0; i<3; ++i)  
         for (size_t j=0; j<3; ++j) {  
           if (eSet.isSecondOrderOperator(i,j)) {  
             generatesElementaryMatrix(PDEOperator::secondorderop, J,  
                                       eSet.secondOrderOperator(i,j),i,j);  
           }  
         }  
     }  
     if (eSet.isDivMuGrad()) {  
       generatesElementaryMatrix(PDEOperator::divmugrad, J, eSet.divMuGrad());  
     }  
   }  
   
   /**  
    * Generates an elementary matrix for a given operator in an  
    * element. The row and column number can be specified when operator  
    * is not scalar: \f$ \partial x_i(w_l)\partial x_j(w_k)\f$ for instance.  
    *  
    * @param operatorType type of the operator  
    * @param J jacobian of the transformation  
    * @param matelem generated elementary matrix  
    * @param i row number  
    * @param j column number  
    */  
   void  
   generatesElementaryMatrix(const PDEOperator::Type operatorType,  
                             const JacobianTransformation& J,  
                             ElementaryMatrixType& matelem,  
                             const size_t i = 0, const size_t j = 0) const  
   {  
     matelem = 0;  
     switch(operatorType) {  
     case PDEOperator::firstorderop: {  
       FiniteElementType::instance().integrateDWjWi(matelem,i,J);  
       matelem *= J.jacobianDet();  
       break;  
     }  
     case PDEOperator::firstorderopTransposed: {  
       FiniteElementType::instance().integrateWjDWi(matelem,i,J);  
       matelem *= J.jacobianDet();  
       break;  
     }  
     case PDEOperator::divmugrad: {  
       FiniteElementType::instance().integrateDWjDWi(matelem,0,0,J);  
       FiniteElementType::instance().integrateDWjDWi(matelem,1,1,J);  
       FiniteElementType::instance().integrateDWjDWi(matelem,2,2,J);  
       matelem *= J.jacobianDet();  
       break;  
     }  
     case PDEOperator::secondorderop: {  
       FiniteElementType::instance().integrateDWjDWi(matelem,i,j,J);  
       matelem *= J.jacobianDet();  
       break;  
     }  
     case PDEOperator::massop: {  
       FiniteElementType::instance().integrateWjWi(matelem,J);  
       matelem *= J.jacobianDet();  
       break;  
     }  
     default: {  
       fferr(2) << '\n' << __FILE__ << ':' << __LINE__  
                << ':' << "Not implemented\n";  
       std::exit(1);  
     }  
     }  
   }  
   
87  public:  public:
88    /**    /**
89     * Assembles the matrix associated to the PDE operators of the PDE     * Assembles the matrix associated to the PDE operators of the PDE
# Line 457  public: Line 314  public:
314      }      }
315    }    }
316    
   /**  
    * Access function to the discretization mesh  
    *  
    * @return the mesh  
    */  
   Mesh& mesh()  
   {  
     return (__mesh);  
   }  
   
   /**  
    * Read only access to the discretization mesh  
    *  
    * @return the mesh  
    */  
   const Mesh& mesh() const  
   {  
     return (__mesh);  
   }  
   
   
317  public:  public:
318    
319    /**    /**
# Line 495  public: Line 331  public:
331                      BaseMatrix& a,                      BaseMatrix& a,
332                      BaseVector& bb,                      BaseVector& bb,
333                      const DegreeOfFreedomSet& dof)                      const DegreeOfFreedomSet& dof)
334      : Discretization(Discretization::FEM,p,a,bb),      : BaseFEMDiscretization<MeshType>(p, m, a, bb, dof)
       __mesh(m),  
       __eSet(problem()),  
       __discretizedOperators(__eSet,problem()),  
       __degreeOfFreedomSet(dof)  
335    {    {
336      ;      ;
337    }    }
# Line 508  public: Line 340  public:
340     * Virtual destructor     * Virtual destructor
341     *     *
342     */     */
343    virtual ~FEMDiscretization()    ~FEMDiscretization()
344    {    {
345      ;      ;
346    }    }
# Line 522  public: Line 354  public:
354   */   */
355  template <>  template <>
356  class FEMDiscretization<Structured3DMesh>  class FEMDiscretization<Structured3DMesh>
357    : public Discretization    : public BaseFEMDiscretization<Structured3DMesh>
358  {  {
359  public:  public:
360    /// The type of mesh used for discretization    /// The type of mesh used for discretization
# Line 549  public: Line 381  public:
381    /// Associated jacobian    /// Associated jacobian
382    typedef FiniteElement::JacobianTransformation   JacobianTransformation;    typedef FiniteElement::JacobianTransformation   JacobianTransformation;
383    
 private:  
   /// Mesh used to perform discretization  
   MeshType& __mesh;  
   
   /// Set of elementary matrices  
   mutable ElementaryMatrixSet <ElementaryMatrixType> __eSet;  
   
   /// Operators that are discretized  
   mutable DiscretizedOperators<ElementaryMatrixType> __discretizedOperators;  
   
   /// Set of degrees of freedom  
   const DegreeOfFreedomSet& __degreeOfFreedomSet;  
   
   /**  
    * Generates elementary vector  
    *  
    * @param eVector the generated elementary vector  
    * @param J the jacobian of the transformation  
    * @param f the function to discretize  
    */  
   void  
   generatesElementaryVector(ElementaryVectorType& eVector,  
                             const JacobianTransformation& J,  
                             const ElementaryVectorType& f) const  
   {  
     FiniteElementType::instance().integrateWj(eVector,J,f);  
   }  
   
   /**  
    * Generates elementary matrices set for a given element  
    *  
    * @param eSet the set of elementary matrices  
    * @param J the jacobian of the transformation  
    */  
   void  
   generatesElementaryMatrix(ElementaryMatrixSet<TinyMatrix<8,8> >& eSet,  
                             const JacobianTransformation& J) const  
   {  
     if (eSet.isMassOperator()) {  
       generatesElementaryMatrix(PDEOperator::massop,  
                                 J, eSet.massOperator());  
     }  
     if (eSet.isFirstOrderOperator()) {  
       for (size_t i=0; i<3; ++i) {  
         if (eSet.isFirstOrderUdxV(i)) {  
           generatesElementaryMatrix(PDEOperator::firstorderopTransposed,  
                                     J,eSet.firstOrderOperatorUdxV(i),i);  
         }  
         if (eSet.isFirstOrderDxUV(i)) {  
           generatesElementaryMatrix(PDEOperator::firstorderop, J,  
                                     eSet.firstOrderOperatorDxUV(i),i);  
         }  
       }  
     }  
     if (eSet.isSecondOrderOperator()) {  
       for (size_t i=0; i<3; ++i)  
         for (size_t j=0; j<3; ++j) {  
           if (eSet.isSecondOrderOperator(i,j)) {  
             generatesElementaryMatrix(PDEOperator::secondorderop, J,  
                                       eSet.secondOrderOperator(i,j),i,j);  
           }  
         }  
     }  
     if (eSet.isDivMuGrad()) {  
       generatesElementaryMatrix(PDEOperator::divmugrad, J, eSet.divMuGrad());  
     }  
   }  
   
   /**  
    * Generates an elementary matrix for a given operator in an  
    * element. The row and column number can be specified when operator  
    * is not scalar: \f$ \partial x_i(w_l)\partial x_j(w_k)\f$ for instance.  
    *  
    * @param operatorType type of the operator  
    * @param J jacobian of the transformation  
    * @param matelem generated elementary matrix  
    * @param i row number  
    * @param j column number  
    */  
   void  
   generatesElementaryMatrix(const PDEOperator::Type operatorType,  
                             const JacobianTransformation& J,  
                             ElementaryMatrixType& matelem,  
                             const size_t i = 0, const size_t j = 0) const  
   {  
     matelem = 0;  
     switch(operatorType) {  
     case PDEOperator::firstorderop: {  
       FiniteElementType::instance().integrateDWjWi(matelem,i,J);  
       matelem *= J.jacobianDet();  
       break;  
     }  
     case PDEOperator::firstorderopTransposed: {  
       FiniteElementType::instance().integrateWjDWi(matelem,i,J);  
       matelem *= J.jacobianDet();  
       break;  
     }  
     case PDEOperator::divmugrad: {  
       FiniteElementType::instance().integrateDWjDWi(matelem,0,0,J);  
       FiniteElementType::instance().integrateDWjDWi(matelem,1,1,J);  
       FiniteElementType::instance().integrateDWjDWi(matelem,2,2,J);  
       matelem *= J.jacobianDet();  
       break;  
     }  
     case PDEOperator::secondorderop: {  
       FiniteElementType::instance().integrateDWjDWi(matelem,i,j,J);  
       matelem *= J.jacobianDet();  
       break;  
     }  
     case PDEOperator::massop: {  
       FiniteElementType::instance().integrateWjWi(matelem,J);  
       matelem *= J.jacobianDet();  
       break;  
     }  
     default: {  
       fferr(2) << '\n' << __FILE__ << ':' << __LINE__  
                << ':' << "Not implemented\n";  
       std::exit(1);  
     }  
     }  
   }  
   
384  public:  public:
385    /**    /**
386     * Assembles the matrix associated to the PDE operators of the PDE     * Assembles the matrix associated to the PDE operators of the PDE
# Line 1150  public: Line 860  public:
860        }        }
861    }    }
862    
   /**  
    * Access function to the discretization mesh  
    *  
    * @return the mesh  
    */  
   Mesh& mesh()  
   {  
     return (__mesh);  
   }  
   
   /**  
    * Read only access to the discretization mesh  
    *  
    * @return the mesh  
    */  
   const Mesh& mesh() const  
   {  
     return (__mesh);  
   }  
   
   
863  public:  public:
864    
865    /**    /**
# Line 1188  public: Line 877  public:
877                      BaseMatrix& a,                      BaseMatrix& a,
878                      BaseVector& bb,                      BaseVector& bb,
879                      const DegreeOfFreedomSet& dof)                      const DegreeOfFreedomSet& dof)
880      : Discretization(Discretization::FEM,p,a,bb),      : BaseFEMDiscretization<Structured3DMesh>(p, m, a, bb, dof)
       __mesh(m),  
       __eSet(problem()),  
       __discretizedOperators(__eSet,problem()),  
       __degreeOfFreedomSet(dof)  
881    {    {
882      ;      ;
883    }    }
884    
885    /**    /**
886     * Virtual destructor     * destructor
887     *     *
888     */     */
889    virtual ~FEMDiscretization()    ~FEMDiscretization()
890    {    {
891      ;      ;
892    }    }

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26