/[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.3 by delpinux, Sat Apr 12 16:16:47 2003 UTC revision 1.4 by delpinux, Sat Apr 12 17:09:16 2003 UTC
# Line 18  Line 18 
18  //  $Id$  //  $Id$
19    
20    
21  #ifndef _FEM_DISCRETIZATION_HPP_  #ifndef FEM_DISCRETIZATION_HPP
22  #define _FEM_DISCRETIZATION_HPP_  #define FEM_DISCRETIZATION_HPP
23    
24    /**
25     * @file   FEMDiscretization.hpp
26     * @author Stephane Del Pino
27     * @date   Sat Apr 12 18:17:18 2003
28     *
29     * @brief  PDE discretization using finite element method.
30     *
31     * Partial Differential Equation discretization using finite element
32     * method. The trick here is that elementary matrices assembling is
33     * factorized in order to optimize the process.
34     *
35     * @warning variational formula is not complete for right hand side
36     * containing partial differential operators.
37     * @todo re-implement rhs discretization using the same design than
38     * for left part.
39     */
40    
41  #include <ElementaryMatrixSet.hpp>  #include <ElementaryMatrixSet.hpp>
42    
# Line 60  Line 76 
76  #include <MeshExpression.hpp>  #include <MeshExpression.hpp>
77    
78    
79  template <typename MeshType>  template <typename GivenMeshType>
80  class FEMDiscretization  class FEMDiscretization
81    : public Discretization    : public Discretization
82  {  {
83  public:  public:
84      /// The type of mesh used for discretization
85      typedef GivenMeshType MeshType;
86    
87      /// The geometry of finite elements
88    typedef typename MeshType::ElementsGeometry ElementsGeometry;    typedef typename MeshType::ElementsGeometry ElementsGeometry;
89    
90    typedef FiniteElementTraits<ElementsGeometry> FiniteElement; /**< The Q1 finite element */    /// The finite element
91      typedef FiniteElementTraits<ElementsGeometry> FiniteElement;
92    
93      /// Finite element type
94    typedef typename FiniteElement::Type FiniteElementType;    typedef typename FiniteElement::Type FiniteElementType;
95    
96      /// Elementary matrices type
97    typedef typename FiniteElement::ElementaryMatrix         ElementaryMatrixType;    typedef typename FiniteElement::ElementaryMatrix         ElementaryMatrixType;
98    
99      /// Elementary vector type
100    typedef typename FiniteElement::ElementaryVector         ElementaryVectorType;    typedef typename FiniteElement::ElementaryVector         ElementaryVectorType;
101    
102      /// type of transformation from reference element
103    typedef typename FiniteElement::Transformation           ConformTransformation;    typedef typename FiniteElement::Transformation           ConformTransformation;
104    
105      /// Associated jacobian
106    typedef typename FiniteElement::JacobianTransformation   JacobianTransformation;    typedef typename FiniteElement::JacobianTransformation   JacobianTransformation;
107    
108  private:  private:
109    //! FEM discretization needs a mesh.    /// Mesh used to perform discretization
110    MeshType& __mesh;    MeshType& __mesh;
111    
112      /// Set of elementary matrices
113    mutable ElementaryMatrixSet <ElementaryMatrixType> __eSet;    mutable ElementaryMatrixSet <ElementaryMatrixType> __eSet;
114    
115      /// Operators that are discretized
116    mutable DiscretizedOperators<ElementaryMatrixType> __discretizedOperators;    mutable DiscretizedOperators<ElementaryMatrixType> __discretizedOperators;
117    
118      /// Set of degrees of freedom
119    const DegreeOfFreedomSet& __degreeOfFreedomSet;    const DegreeOfFreedomSet& __degreeOfFreedomSet;
120    
121      /**
122       * Generates elementary vector
123       *
124       * @param eVector the generated elementary vector
125       * @param J the jacobian of the transformation
126       * @param f the function to discretize
127       */
128    void    void
129    generatesElementaryVector(ElementaryVectorType& eVector,    generatesElementaryVector(ElementaryVectorType& eVector,
130                              const JacobianTransformation& J,                              const JacobianTransformation& J,
# Line 94  private: Line 133  private:
133      FiniteElementType::instance().integrateWj(eVector,J,f);      FiniteElementType::instance().integrateWj(eVector,J,f);
134    }    }
135    
136      /**
137       * Generates elementary matrices set for a given element
138       *
139       * @param eSet the set of elementary matrices
140       * @param J the jacobian of the transformation
141       */
142    void    void
143    generatesElementaryMatrix(ElementaryMatrixSet<ElementaryMatrixType>& eSet,    generatesElementaryMatrix(ElementaryMatrixSet<ElementaryMatrixType>& eSet,
144                              const JacobianTransformation& J) const                              const JacobianTransformation& J) const
# Line 128  private: Line 173  private:
173      }      }
174    }    }
175    
176      /**
177       * Generates an elementary matrix for a given operator in an
178       * element. The row and column number can be specified when operator
179       * is not scalar: \f$ \partial x_i(w_l)\partial x_j(w_k)\f$ for instance.
180       *
181       * @param operatorType type of the operator
182       * @param J jacobian of the transformation
183       * @param matelem generated elementary matrix
184       * @param i row number
185       * @param j column number
186       */
187    void    void
188    generatesElementaryMatrix(const PDEOperator::Type operatorType,    generatesElementaryMatrix(const PDEOperator::Type operatorType,
189                              const JacobianTransformation& J,                              const JacobianTransformation& J,
# Line 172  private: Line 227  private:
227      }      }
228    }    }
229    
   
230  public:  public:
231    /*! Assembles the matrix associated to the PDE operators of the PDE    /**
232      problem.     * Assembles the matrix associated to the PDE operators of the PDE
233       * problem.
234       *
235     */     */
236    void assembleMatrix()    void assembleMatrix()
237    {    {
# Line 243  public: Line 299  public:
299      ffout(3) << "Assembling cost: " << t << '\n';      ffout(3) << "Assembling cost: " << t << '\n';
300    }    }
301    
302    /*! Applies directly the operator discretization to the vector X.    /**
303       *  Applies directly the operator discretization to the vector X.
304       *
305       * @param X input vector
306       * @param Z \f$ z=Ax \f$
307     */     */
308    void timesX(const BaseVector& X, BaseVector& Z) const    void timesX(const BaseVector& X, BaseVector& Z) const
309    {    {
# Line 290  public: Line 350  public:
350      }      }
351    }    }
352    
353    /*! Gets directly the diagonal of the discretization to the vector X.    /**
354       * Computes diagonal of the operator
355       *
356       * @param Z diagonal of the operator
357     */     */
358    void getDiagonal(BaseVector& Z) const    void getDiagonal(BaseVector& Z) const
359    {    {
# Line 333  public: Line 396  public:
396      }      }
397    }    }
398    
399    /*! Assembles the matrix associated to the PDE operators of the PDE    /**
400      problem.     * Second member assembling
401       *
402     */     */
403    void assembleSecondMember()    void assembleSecondMember()
404    {    {
# Line 393  public: Line 457  public:
457      }      }
458    }    }
459    
460      /**
461    //! mesh     * Access function to the discretization mesh
462       *
463       * @return the mesh
464       */
465    Mesh& mesh()    Mesh& mesh()
466    {    {
467      return (__mesh);      return (__mesh);
468    }    }
469    
470    //! mesh    /**
471       * Read only access to the discretization mesh
472       *
473       * @return the mesh
474       */
475    const Mesh& mesh() const    const Mesh& mesh() const
476    {    {
477      return (__mesh);      return (__mesh);
# Line 409  public: Line 480  public:
480    
481  public:  public:
482    
483    //! Constructor    /**
484       * Constructor of the discretization
485       *
486       * @param p the problem
487       * @param m the mesh used for discretization
488       * @param a matrix storing discretization
489       * @param bb vector that stores second member discretization
490       * @param dof degrees of freedom set
491       *
492       */
493    FEMDiscretization(const Problem& p,    FEMDiscretization(const Problem& p,
494                      MeshType& m,                      MeshType& m,
495                      BaseMatrix& a,                      BaseMatrix& a,
# Line 424  public: Line 504  public:
504      ;      ;
505    }    }
506    
507    //! Virtual destructor.    /**
508       * Virtual destructor
509       *
510       */
511    virtual ~FEMDiscretization()    virtual ~FEMDiscretization()
512    {    {
513      ;      ;
# Line 432  public: Line 515  public:
515  };  };
516    
517    
518    /**
519     * This is specialisation in the particular case of Cartesian grids.
520     * It is an important specialization since the elementary matrices are
521     * computed \bf once for all!
522     */
523  template <>  template <>
524  class FEMDiscretization<Structured3DMesh>  class FEMDiscretization<Structured3DMesh>
525    : public Discretization    : public Discretization
526  {  {
527  public:  public:
528      /// The type of mesh used for discretization
529    typedef Structured3DMesh MeshType;    typedef Structured3DMesh MeshType;
530    
531      /// The geometry of finite elements
532    typedef MeshType::ElementsGeometry ElementsGeometry;    typedef MeshType::ElementsGeometry ElementsGeometry;
533    
534      /// The finite element
535    typedef FiniteElementTraits<ElementsGeometry> FiniteElement; /**< The Q1 finite element */    typedef FiniteElementTraits<ElementsGeometry> FiniteElement; /**< The Q1 finite element */
536    
537      /// Finite element type
538    typedef FiniteElement::Type FiniteElementType;    typedef FiniteElement::Type FiniteElementType;
539    
540      /// Elementary matrices type
541    typedef FiniteElement::ElementaryMatrix         ElementaryMatrixType;    typedef FiniteElement::ElementaryMatrix         ElementaryMatrixType;
542    
543      /// Elementary vector type
544    typedef FiniteElement::ElementaryVector         ElementaryVectorType;    typedef FiniteElement::ElementaryVector         ElementaryVectorType;
545    
546      /// type of transformation from reference element
547    typedef FiniteElement::Transformation           ConformTransformation;    typedef FiniteElement::Transformation           ConformTransformation;
548    
549      /// Associated jacobian
550    typedef FiniteElement::JacobianTransformation   JacobianTransformation;    typedef FiniteElement::JacobianTransformation   JacobianTransformation;
551    
552  private:  private:
553    //! FEM discretization needs a mesh.    /// Mesh used to perform discretization
554    MeshType& __mesh;    MeshType& __mesh;
555    
556      /// Set of elementary matrices
557    mutable ElementaryMatrixSet <ElementaryMatrixType> __eSet;    mutable ElementaryMatrixSet <ElementaryMatrixType> __eSet;
558    
559      /// Operators that are discretized
560    mutable DiscretizedOperators<ElementaryMatrixType> __discretizedOperators;    mutable DiscretizedOperators<ElementaryMatrixType> __discretizedOperators;
561    
562      /// Set of degrees of freedom
563    const DegreeOfFreedomSet& __degreeOfFreedomSet;    const DegreeOfFreedomSet& __degreeOfFreedomSet;
564    
565      /**
566       * Generates elementary vector
567       *
568       * @param eVector the generated elementary vector
569       * @param J the jacobian of the transformation
570       * @param f the function to discretize
571       */
572    void    void
573    generatesElementaryVector(ElementaryVectorType& eVector,    generatesElementaryVector(ElementaryVectorType& eVector,
574                              const JacobianTransformation& J,                              const JacobianTransformation& J,
# Line 469  private: Line 577  private:
577      FiniteElementType::instance().integrateWj(eVector,J,f);      FiniteElementType::instance().integrateWj(eVector,J,f);
578    }    }
579    
580      /**
581       * Generates elementary matrices set for a given element
582       *
583       * @param eSet the set of elementary matrices
584       * @param J the jacobian of the transformation
585       */
586    void    void
587    generatesElementaryMatrix(ElementaryMatrixSet<TinyMatrix<8,8> >& eSet,    generatesElementaryMatrix(ElementaryMatrixSet<TinyMatrix<8,8> >& eSet,
588                              const JacobianTransformation& J) const                              const JacobianTransformation& J) const
# Line 503  private: Line 617  private:
617      }      }
618    }    }
619    
620      /**
621       * Generates an elementary matrix for a given operator in an
622       * element. The row and column number can be specified when operator
623       * is not scalar: \f$ \partial x_i(w_l)\partial x_j(w_k)\f$ for instance.
624       *
625       * @param operatorType type of the operator
626       * @param J jacobian of the transformation
627       * @param matelem generated elementary matrix
628       * @param i row number
629       * @param j column number
630       */
631    void    void
632    generatesElementaryMatrix(const PDEOperator::Type operatorType,    generatesElementaryMatrix(const PDEOperator::Type operatorType,
633                              const JacobianTransformation& J,                              const JacobianTransformation& J,
# Line 546  private: Line 671  private:
671      }      }
672    }    }
673    
   
674  public:  public:
675    /*! Assembles the matrix associated to the PDE operators of the PDE    /**
676      problem.     * Assembles the matrix associated to the PDE operators of the PDE
677       * problem.
678       *
679     */     */
680    void assembleMatrix()    void assembleMatrix()
681    {    {
# Line 618  public: Line 744  public:
744      ffout(3) << "Assembling cost: " << t << '\n';      ffout(3) << "Assembling cost: " << t << '\n';
745    }    }
746    
747    /*! Applies directly the operator discretization to the vector X.    /**
748       *  Applies directly the operator discretization to the vector X.
749       *
750       * @param X input vector
751       * @param Z \f$ z=Ax \f$
752     */     */
753    void timesX(const BaseVector& X, BaseVector& Z) const    void timesX(const BaseVector& X, BaseVector& Z) const
754    {    {
# Line 665  public: Line 795  public:
795      }      }
796    }    }
797    
798    /*! Gets directly the diagonal of the discretization to the vector X.    /**
799       * Computes diagonal of the operator
800       *
801       * @param Z diagonal of the operator
802     */     */
803    void getDiagonal(BaseVector& Z) const    void getDiagonal(BaseVector& Z) const
804    {    {
# Line 706  public: Line 839  public:
839      }      }
840    }    }
841    
842    /*! Assembles the matrix associated to the PDE operators of the PDE    /**
843      problem.     * Second member assembling
844       *
845     */     */
846    void assembleSecondMember()    void assembleSecondMember()
847    {    {
# Line 1016  public: Line 1150  public:
1150        }        }
1151    }    }
1152    
1153    //! mesh    /**
1154       * Access function to the discretization mesh
1155       *
1156       * @return the mesh
1157       */
1158    Mesh& mesh()    Mesh& mesh()
1159    {    {
1160      return (__mesh);      return (__mesh);
1161    }    }
1162    
1163    //! mesh    /**
1164       * Read only access to the discretization mesh
1165       *
1166       * @return the mesh
1167       */
1168    const Mesh& mesh() const    const Mesh& mesh() const
1169    {    {
1170      return (__mesh);      return (__mesh);
# Line 1031  public: Line 1173  public:
1173    
1174  public:  public:
1175    
1176    //! Constructor    /**
1177       * Constructor of the discretization
1178       *
1179       * @param p the problem
1180       * @param m the mesh used for discretization
1181       * @param a matrix storing discretization
1182       * @param bb vector that stores second member discretization
1183       * @param dof degrees of freedom set
1184       *
1185       */
1186    FEMDiscretization(const Problem& p,    FEMDiscretization(const Problem& p,
1187                      Structured3DMesh& m,                      Structured3DMesh& m,
1188                      BaseMatrix& a,                      BaseMatrix& a,
# Line 1046  public: Line 1197  public:
1197      ;      ;
1198    }    }
1199    
1200    //! Virtual destructor.    /**
1201       * Virtual destructor
1202       *
1203       */
1204    virtual ~FEMDiscretization()    virtual ~FEMDiscretization()
1205    {    {
1206      ;      ;
# Line 1054  public: Line 1208  public:
1208  };  };
1209    
1210    
1211  #endif // _FEM_DISCRETIZATION_HPP_  #endif // FEM_DISCRETIZATION_HPP

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

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