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 |
|
|
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, |
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 |
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, |
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 |
{ |
{ |
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 |
{ |
{ |
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 |
{ |
{ |
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 |
{ |
{ |
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); |
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, |
504 |
; |
; |
505 |
} |
} |
506 |
|
|
507 |
//! Virtual destructor. |
/** |
508 |
|
* Virtual destructor |
509 |
|
* |
510 |
|
*/ |
511 |
virtual ~FEMDiscretization() |
virtual ~FEMDiscretization() |
512 |
{ |
{ |
513 |
; |
; |
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, |
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 |
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, |
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 |
{ |
{ |
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 |
{ |
{ |
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 |
{ |
{ |
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 |
{ |
{ |
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); |
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, |
1197 |
; |
; |
1198 |
} |
} |
1199 |
|
|
1200 |
//! Virtual destructor. |
/** |
1201 |
|
* Virtual destructor |
1202 |
|
* |
1203 |
|
*/ |
1204 |
virtual ~FEMDiscretization() |
virtual ~FEMDiscretization() |
1205 |
{ |
{ |
1206 |
; |
; |
1208 |
}; |
}; |
1209 |
|
|
1210 |
|
|
1211 |
#endif // _FEM_DISCRETIZATION_HPP_ |
#endif // FEM_DISCRETIZATION_HPP |