FEMDiscretization< Structured3DMesh, TypeOfDiscretization > Class Template Reference

#include <FEMDiscretization.hpp>

Inheritance diagram for FEMDiscretization< Structured3DMesh, TypeOfDiscretization >:

Inheritance graph
[legend]
Collaboration diagram for FEMDiscretization< Structured3DMesh, TypeOfDiscretization >:

Collaboration graph
[legend]

List of all members.

Public Types

typedef Structured3DMesh MeshType
 The type of mesh used for discretization.
typedef MeshType::CellType CellType
 The geometry of finite elements.
typedef FiniteElementTraits
< CellType,
TypeOfDiscretization >::Type 		FiniteElement
 The finite element.
typedef
FiniteElement::QuadratureType 		QuadratureType
 quadrature type
typedef
FiniteElement::ElementaryMatrix 		ElementaryMatrixType
 Elementary matrices type.
typedef
FiniteElement::ElementaryVector 		ElementaryVectorType
 Elementary vector type.
typedef FiniteElementTraits
< CellType,
TypeOfDiscretization >
::Transformation 		ConformTransformation
 type of transformation from reference element
typedef FiniteElementTraits
< CellType,
TypeOfDiscretization >
::JacobianTransformation 		JacobianTransformation
 Associated jacobian.

Public Member Functions

void assembleMatrix ()
void timesX (const BaseVector &X, BaseVector &Z) const
void transposedTimesX (const BaseVector &X, BaseVector &Z) const
void getDiagonal (BaseVector &Z) const
void assembleSecondMember ()
 FEMDiscretization (const Problem &p, const Structured3DMesh &m, const DiscretizationType &discretisationType, BaseMatrix &a, BaseVector &bb, const DegreeOfFreedomSet &dof)
 ~FEMDiscretization ()
void setDirichletList (const Vector< bool > &dirichletList)
const Problem & problem () const
BaseMatrix & A ()
BaseVector & b ()

Protected Member Functions

void generatesElementaryVector (ElementaryVectorType &eVector, const JacobianTransformation &J, const TinyVector< FiniteElement::numberOfQuadraturePoints, real_t > &f) const
void generatesElementaryMatrix (ElementaryMatrixSet< ElementaryMatrixType > &eSet, const JacobianTransformation &J) const
void generatesElementaryMatrix (const PDEOperator::Type operatorType, const JacobianTransformation &J, ElementaryMatrixType &matelem, const size_t i=0, const size_t j=0) const

Protected Attributes

const MeshType & __mesh
 Mesh used to perform discretization.
ElementaryMatrixSet
< ElementaryMatrixType > 		__eSet
 Set of elementary matrices.
DiscretizedOperators
< ElementaryMatrixType > 		__discretizedOperators
 Operators that are discretized.
const DegreeOfFreedomSet & __degreeOfFreedomSet
 Set of degrees of freedom.
const DiscretizationType __discretizationType
const Problem & __problem
 The PDEProblem to discretize.
BaseMatrix & __A
 The matrix which will contain the discretization.
BaseVector & __b
 The second member.
const Vector< bool > * __dirichletList
 elimination dirichlet informations

Detailed Description

template<>
class FEMDiscretization< Structured3DMesh, TypeOfDiscretization >

This is specialisation in the particular case of Cartesian grids. It is an important specialization since the elementary matrices are computed once for all!

Definition at line 654 of file FEMDiscretization.hpp.

Member Typedef Documentation

typedef Structured3DMesh FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::MeshType

The type of mesh used for discretization.

Reimplemented from BaseFEMDiscretization< Structured3DMesh, TypeOfDiscretization >.

Definition at line 660 of file FEMDiscretization.hpp.

typedef MeshType::CellType FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::CellType

The geometry of finite elements.

Reimplemented from BaseFEMDiscretization< Structured3DMesh, TypeOfDiscretization >.

Definition at line 663 of file FEMDiscretization.hpp.

typedef FiniteElementTraits<CellType, TypeOfDiscretization>::Type FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::FiniteElement

The finite element.

Reimplemented from BaseFEMDiscretization< Structured3DMesh, TypeOfDiscretization >.

Definition at line 667 of file FEMDiscretization.hpp.

typedef FiniteElement::QuadratureType FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::QuadratureType

quadrature type

Definition at line 670 of file FEMDiscretization.hpp.

typedef FiniteElement::ElementaryMatrix FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::ElementaryMatrixType

Elementary matrices type.

Reimplemented from BaseFEMDiscretization< Structured3DMesh, TypeOfDiscretization >.

Definition at line 673 of file FEMDiscretization.hpp.

typedef FiniteElement::ElementaryVector FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::ElementaryVectorType

Elementary vector type.

Reimplemented from BaseFEMDiscretization< Structured3DMesh, TypeOfDiscretization >.

Definition at line 676 of file FEMDiscretization.hpp.

typedef FiniteElementTraits<CellType, TypeOfDiscretization>::Transformation FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::ConformTransformation

type of transformation from reference element

Reimplemented from BaseFEMDiscretization< Structured3DMesh, TypeOfDiscretization >.

Definition at line 680 of file FEMDiscretization.hpp.

typedef FiniteElementTraits<CellType, TypeOfDiscretization>::JacobianTransformation FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::JacobianTransformation

Associated jacobian.

Reimplemented from BaseFEMDiscretization< Structured3DMesh, TypeOfDiscretization >.

Definition at line 684 of file FEMDiscretization.hpp.

Constructor & Destructor Documentation

FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::FEMDiscretization ( const Problem &  p,
const Structured3DMesh &  m,
const DiscretizationType &  discretisationType,
BaseMatrix &  a,
BaseVector &  bb,
const DegreeOfFreedomSet &  dof 
) [inline]

Constructor of the discretization

Parameters:
p the problem
m the mesh used for discretization
discretizationType the type of discretization
a matrix storing discretization
bb vector that stores second member discretization
dof degrees of freedom set

Definition at line 1298 of file FEMDiscretization.hpp.

References Discretization::__discretizationType, StaticBase< SolverInformationCenter >::instance(), SolverInformationCenter::pushDiscretizationType(), and SolverInformationCenter::pushMesh().

01304     : BaseFEMDiscretization<Structured3DMesh,
01305                             TypeOfDiscretization>(p, m, discretisationType, a, bb, dof)
01306   {
01307     SolverInformationCenter::instance().pushMesh(&m);
01308     SolverInformationCenter::instance().pushDiscretizationType(&(this->__discretizationType));
01309   }

Here is the call graph for this function:

FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::~FEMDiscretization (  )  [inline]

destructor

Definition at line 1315 of file FEMDiscretization.hpp.

References StaticBase< SolverInformationCenter >::instance(), and SolverInformationCenter::pop().

01316   {
01317     SolverInformationCenter::instance().pop();
01318   }

Here is the call graph for this function:

Member Function Documentation

void FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::assembleMatrix (  )  [inline, virtual]

Assembles the matrix associated to the PDE operators of the PDE problem.

Implements Discretization.

Definition at line 692 of file FEMDiscretization.hpp.

References Discretization::__A, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__degreeOfFreedomSet, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__discretizedOperators, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__eSet, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__mesh, Discretization::A(), BaseMatrix::doubleHashedMatrix, ffout(), BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::generatesElementaryMatrix(), DegreeOfFreedomSet::positionsSet(), UnAssembledMatrix::setDiscretization(), Timer::start(), Timer::stop(), BaseMatrix::type(), BaseMatrix::unAssembled, and ErrorHandler::unexpected.

00693   {
00694     const ScalarDegreeOfFreedomPositionsSet& dofPositions
00695       = this->__degreeOfFreedomSet.positionsSet(0);
00696     
00697     switch (this->__A.type()) {
00698     case BaseMatrix::doubleHashedMatrix: {
00699       Timer t;
00700       t.start();
00701 
00702       ffout(2) << "- assembling matrix\n";
00703 
00704       DoubleHashedMatrix& A = dynamic_cast<DoubleHashedMatrix&>(this->__A);
00705 
00706       MeshType::const_iterator icell0(this->__mesh);
00707       ConformTransformation T0(*icell0);
00708       JacobianTransformation J(T0);
00709       generatesElementaryMatrix(this->__eSet,J);
00710 
00711       for(MeshType::const_iterator icell(this->__mesh);
00712           not(icell.end()); ++icell) {
00713         const CellType& K = *icell;
00714         ConformTransformation T(K);
00715         TinyVector<3> massCenter;
00716         T.value(FiniteElement::massCenter(), massCenter);
00717 
00718         TinyVector<FiniteElement::numberOfDegreesOfFreedom,
00719           size_t> index;
00720         const size_t cellNumber = icell.number();
00721 
00722         for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00723           index[l] = dofPositions(cellNumber,l);
00724         }
00725         for(typename DiscretizedOperators<ElementaryMatrixType>
00726               ::iterator iOperator = this->__discretizedOperators.begin();
00727             iOperator != this->__discretizedOperators.end(); ++iOperator) {
00728 
00729           const real_t coef = ((*iOperator).second)(massCenter);
00730 
00731           const ElementaryMatrixType& elementaryMatrix
00732             = (*(*iOperator).first);
00733 
00734           for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
00735             const size_t I1 // Index of Ith value
00736               = this->__degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
00737             for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
00738               const size_t I2 // Index of Jth value
00739                 = this->__degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
00740 
00741               A(I1,I2) += coef*elementaryMatrix(l,m);
00742             }
00743           }
00744         }
00745       }
00746       ffout(2) << "- assembling matrix: done";
00747       t.stop();
00748       ffout(3) << " [cost: " << t << ']';
00749       ffout(2) << '\n';
00750 
00751       break;
00752     }
00753     case BaseMatrix::unAssembled: {
00754       UnAssembledMatrix& A = dynamic_cast<UnAssembledMatrix&>(this->__A);
00755       A.setDiscretization(this);
00756       break;
00757     }
00758 
00759     default: {
00760       throw ErrorHandler(__FILE__,__LINE__,
00761                          "unexpected matrix type",
00762                          ErrorHandler::unexpected);
00763     }
00764     }
00765   }

Here is the call graph for this function:

void FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::timesX ( const BaseVector &  X,
BaseVector &  Z 
) const [inline, virtual]

Applies directly the operator discretization to the vector X.

Parameters:
X input vector
Z $ z=Ax $

Implements Discretization.

Definition at line 773 of file FEMDiscretization.hpp.

References BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__degreeOfFreedomSet, Discretization::__dirichletList, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__discretizedOperators, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__eSet, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__mesh, Mesh::T_iterator< MeshType, CellType >::end(), BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::generatesElementaryMatrix(), Mesh::T_iterator< MeshType, CellType >::number(), and DegreeOfFreedomSet::positionsSet().

00774   {
00775     const Vector<real_t>& x = dynamic_cast<const Vector<real_t>&>(X);
00776 
00777     Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
00778     z=0;
00779 
00780     const ScalarDegreeOfFreedomPositionsSet& dofPositions
00781       = this->__degreeOfFreedomSet.positionsSet(0);
00782   
00783     MeshType::const_iterator icell(this->__mesh);
00784     ConformTransformation T0(*icell);
00785     JacobianTransformation J(T0);
00786     generatesElementaryMatrix(this->__eSet,J);
00787 
00788     for(MeshType::const_iterator icell(this->__mesh);
00789         not(icell.end()); ++icell) {
00790       const CellType& K = *icell;
00791       ConformTransformation T(K);
00792 
00793       const size_t cellNumber = icell.number();
00794 
00795       TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
00796       for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00797         index[l] = dofPositions(cellNumber,l);
00798       }
00799 
00800       TinyVector<3> massCenter;
00801       T.value(FiniteElement::massCenter(), massCenter);
00802 
00803       for(typename DiscretizedOperators<ElementaryMatrixType>
00804             ::iterator iOperator = this->__discretizedOperators.begin();
00805           iOperator != this->__discretizedOperators.end(); ++iOperator) {
00806         const real_t coef = ((*iOperator).second)(massCenter);
00807 
00808         const ElementaryMatrixType& elementaryMatrix = (*(*iOperator).first);
00809 
00810 #warning can optimize this loops by precomputing (i and) j
00811         for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
00812           // Number of the degree of freedom
00813           const size_t I1 = this->__degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
00814           if (not((*this->__dirichletList)[I1])) {
00815             for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
00816               // Number of the degree of freedom
00817               const size_t I2 = this->__degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
00818               if (not((*this->__dirichletList)[I2])) {
00819                 z[I1] += coef*elementaryMatrix(l,m)*x[I2];
00820               }
00821             }
00822           }
00823         }
00824       }
00825     }
00826   }

Here is the call graph for this function:

void FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::transposedTimesX ( const BaseVector &  X,
BaseVector &  Z 
) const [inline, virtual]

Applies directly the operator discretization to the vector X.

Parameters:
X input vector
Z $ z=A^T x $

Implements Discretization.

Definition at line 834 of file FEMDiscretization.hpp.

References BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__degreeOfFreedomSet, Discretization::__dirichletList, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__discretizedOperators, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__eSet, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__mesh, Mesh::T_iterator< MeshType, CellType >::end(), BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::generatesElementaryMatrix(), Mesh::T_iterator< MeshType, CellType >::number(), and DegreeOfFreedomSet::positionsSet().

00835   {
00836     const Vector<real_t>& x = dynamic_cast<const Vector<real_t>&>(X);
00837 
00838     Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
00839     z=0;
00840 
00841     const ScalarDegreeOfFreedomPositionsSet& dofPositions
00842       = this->__degreeOfFreedomSet.positionsSet(0);
00843   
00844     MeshType::const_iterator icell(this->__mesh);
00845     ConformTransformation T0(*icell);
00846     JacobianTransformation J(T0);
00847     generatesElementaryMatrix(this->__eSet,J);
00848 
00849     for(MeshType::const_iterator icell(this->__mesh);
00850         not(icell.end()); ++icell) {
00851       const CellType& K = *icell;
00852       ConformTransformation T(K);
00853 
00854       const size_t cellNumber = icell.number();
00855 
00856       TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
00857       for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00858         index[l] = dofPositions(cellNumber,l);
00859       }
00860 
00861       TinyVector<3> massCenter;
00862       T.value(FiniteElement::massCenter(), massCenter);
00863 
00864       for(typename DiscretizedOperators<ElementaryMatrixType>
00865             ::iterator iOperator = this->__discretizedOperators.begin();
00866           iOperator != this->__discretizedOperators.end(); ++iOperator) {
00867         const real_t coef = ((*iOperator).second)(massCenter);
00868 
00869         const ElementaryMatrixType& elementaryMatrix = (*(*iOperator).first);
00870 
00871 #warning can optimize this loops by precomputing (i and) j
00872         for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
00873           // Number of the degree of freedom
00874           const size_t I1 = this->__degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
00875           if (not((*this->__dirichletList)[I1])) {
00876             for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
00877               // Number of the degree of freedom
00878               const size_t I2 = this->__degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
00879               if (not((*this->__dirichletList)[I2])) {
00880                 z[I2] += coef*elementaryMatrix(l,m)*x[I1];
00881               }
00882             }
00883           }
00884         }
00885       }
00886     }
00887   }

Here is the call graph for this function:

void FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::getDiagonal ( BaseVector &  Z  )  const [inline, virtual]

Computes diagonal of the operator

Parameters:
Z diagonal of the operator

Implements Discretization.

Definition at line 894 of file FEMDiscretization.hpp.

References BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__degreeOfFreedomSet, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__discretizedOperators, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__eSet, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__mesh, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::generatesElementaryMatrix(), and DegreeOfFreedomSet::positionsSet().

00895   {
00896     Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
00897 
00898     const ScalarDegreeOfFreedomPositionsSet& dofPositions
00899       = this->__degreeOfFreedomSet.positionsSet(0);
00900   
00901     MeshType::const_iterator icell0(this->__mesh);
00902     ConformTransformation T0(*icell0);
00903     JacobianTransformation J(T0);
00904     generatesElementaryMatrix(this->__eSet,J);
00905 
00906     for(MeshType::const_iterator icell(this->__mesh);
00907         not(icell.end()); ++icell) {
00908       const CellType& K = *icell;
00909       ConformTransformation T(K);
00910 
00911       const size_t cellNumber = icell.number();
00912 
00913       TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
00914       for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00915         index[l] = dofPositions(cellNumber,l);
00916       }
00917 
00918       TinyVector<3> massCenter;
00919       T.value(FiniteElement::massCenter(), massCenter);
00920 
00921       for(typename DiscretizedOperators<ElementaryMatrixType>
00922             ::iterator iOperator = this->__discretizedOperators.begin();
00923           iOperator != this->__discretizedOperators.end(); ++iOperator) {
00924         // Only diagonal opertors have to be considerate
00925         if (((*iOperator).second).i() == ((*iOperator).second).j()) {
00926           const real_t coef = ((*iOperator).second)(massCenter);
00927         
00928           const ElementaryMatrixType& elementaryMatrix = (*(*iOperator).first);
00929           for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
00930             const size_t I = // Index of Ith value
00931               this->__degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
00932             z[I] += coef*elementaryMatrix(l,l);
00933           }
00934         }
00935       }
00936     }
00937   }

Here is the call graph for this function:

void FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::assembleSecondMember (  )  [inline, virtual]

Second member assembling

The elementary vector

Implements Discretization.

Definition at line 943 of file FEMDiscretization.hpp.

References Discretization::__b, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__degreeOfFreedomSet, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__mesh, Discretization::b(), FEMFunctionBase::baseMesh(), VariationalProblem::beginLinearOperator(), Mesh::cartesianHexahedraMesh, PDEOperator::divmugrad, Mesh::T_iterator< MeshType, CellType >::end(), VariationalProblem::endLinearOperator(), VariationalOperatorFgradGgradV::f(), VariationalOperatorFdxV::f(), VariationalOperatorFdxGV::f(), VariationalOperatorFV::f(), VariationalLinearOperator::FdxGV, VariationalLinearOperator::FdxV, ScalarFunctionBase::femfunction, ffout(), VariationalLinearOperator::FgradGgradV, Structured3DMesh::find(), PDEOperator::firstorderop, PDEOperator::firstorderopTransposed, VariationalLinearOperator::FV, VariationalOperatorFgradGgradV::g(), VariationalOperatorFdxGV::g(), BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::generatesElementaryMatrix(), BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::generatesElementaryVector(), FEMFunctionBase::gradient(), VariationalProblem::hasJumpOrMean(), Cell::isFictitious(), ErrorHandler::normal, Mesh::T_iterator< MeshType, CellType >::number(), VariationalOperatorFdxV::number(), VariationalOperatorFdxGV::number(), ScalarDegreeOfFreedomPositionsSet::number(), PDESystem::numberOfUnknown(), Problem::pde, DegreeOfFreedomSet::positionsSet(), Discretization::problem(), PDESystem::secondMember(), Timer::start(), Timer::stop(), VariationalOperator::testFunctionNumber(), ScalarFunctionBase::type(), ErrorHandler::unexpected, Problem::variational, ScalarDegreeOfFreedomPositionsSet::vertex(), and Cell::volume().

00944   {
00945     Timer t;
00946     t.start();
00947 
00948     ffout(2) << "- assembling second member\n";
00949 
00950     const ScalarDegreeOfFreedomPositionsSet& dofPositions
00951       = this->__degreeOfFreedomSet.positionsSet(0);
00952   
00954     ElementaryVectorType eVect;
00955     Vector<real_t>& b = (static_cast<Vector<real_t>&>(this->__b));
00956     b = 0;
00957 
00958     const QuadratureType& referenceQuadrature = QuadratureType::instance();
00959 
00960     switch (this->problem().type()) {
00961     case Problem::pde: {
00962       const PDESystem& pdeSystem
00963         = static_cast<const PDESystem&>(this->problem());
00964 
00965       for(MeshType::const_iterator icell(this->__mesh);
00966             not(icell.end()); ++icell) {
00967         const CellType& K = *icell;
00968         eVect = 0;
00969         ConformTransformation T(K);
00970         JacobianTransformation J(T);
00971         TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
00972         TinyVector<QuadratureType::numberOfQuadraturePoints,
00973                    TinyVector<3> > quadrature;
00974 
00975         for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
00976           T.value(referenceQuadrature[l], quadrature[l]);
00977         }
00978 
00979         const size_t cellNumber = icell.number();
00980         for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00981           index[l] = this->__degreeOfFreedomSet.positionsSet(0)(cellNumber,l);
00982         }
00983 
00984         for (size_t i=0; i<pdeSystem.numberOfUnknown(); ++i) {
00985           const ScalarFunctionBase& F = pdeSystem.secondMember(i);
00986 
00987           TinyVector<QuadratureType::numberOfQuadraturePoints> f;
00988 
00989           for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
00990             f[l] = F(quadrature[l]);
00991           }
00992 
00993           generatesElementaryVector(eVect,J,f);
00994 
00995           for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00996             const size_t dof = this->__degreeOfFreedomSet(i,index[l]);
00997             b[dof] += eVect[l]*J.jacobianDet();
00998           }
00999         }
01000       }
01001       break;
01002     }
01003     case Problem::variational: {
01004       const VariationalProblem& variationalProblem
01005         =  dynamic_cast<const VariationalProblem&>(this->problem());
01006 
01007       if (variationalProblem.hasJumpOrMean()) {
01008         throw ErrorHandler(__FILE__,__LINE__,
01009                            "cannot discretize problems with jumps or means in FEM !",
01010                            ErrorHandler::normal);
01011       }
01012 
01013       for (VariationalProblem::linearOperatorConst_iterator i
01014              = variationalProblem.beginLinearOperator();
01015            i != variationalProblem.endLinearOperator(); ++i) {
01016         switch ((*(*i)).type()) {
01017         case VariationalLinearOperator::FV: {
01018           const VariationalOperatorFV& fv
01019             = dynamic_cast<const VariationalOperatorFV&>(*(*i));
01020 
01021           const ScalarFunctionBase& F = fv.f();
01022           const size_t testFunctionNumber = fv.testFunctionNumber();
01023 
01024           for(MeshType::const_iterator icell(this->__mesh);
01025               not(icell.end()); ++icell) {
01026             const CellType& K = *icell;
01027             eVect = 0;
01028             ConformTransformation T(K);
01029             JacobianTransformation J(T);
01030 
01031             // computes local quadrature vertex
01032             TinyVector<QuadratureType::numberOfQuadraturePoints,
01033                        TinyVector<3, real_t> > quadrature;
01034             for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
01035               T.value(referenceQuadrature[l], quadrature[l]);
01036             }
01037 
01038             const size_t cellNumber = icell.number();
01039 
01040             TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
01041             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
01042               index[l] = dofPositions(cellNumber,l);
01043             }
01044 
01045             ElementaryVectorType f;
01046 
01047             for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
01048               f[l] = F(quadrature[l]);
01049             }
01050 
01051             generatesElementaryVector(eVect,J,f);
01052 
01053             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
01054               const size_t dof = this->__degreeOfFreedomSet(testFunctionNumber,index[l]);
01055               b[dof] += eVect[l]*J.jacobianDet();
01056             }
01057           }
01058           break;
01059         }
01060         case VariationalLinearOperator::FdxGV: {
01061           const VariationalOperatorFdxGV& fv
01062             = dynamic_cast<const VariationalOperatorFdxGV&>(*(*i));
01063 
01064           const ScalarFunctionBase& F = fv.f();
01065           const ScalarFunctionBase& G = fv.g();
01066           const size_t testFunctionNumber = fv.testFunctionNumber();
01067 
01068           Vector<real_t> gValues (dofPositions.number());
01069 
01070           for (size_t i=0; i<dofPositions.number(); ++i) {
01071             gValues[i] = G(dofPositions.vertex(i));
01072           }
01073 
01074           MeshType::const_iterator icell0(this->__mesh);
01075           ConformTransformation T0(*icell0);
01076           JacobianTransformation J(T0);
01077 
01078           ElementaryMatrixType edxUV = 0;
01079           generatesElementaryMatrix(PDEOperator::firstorderop,
01080                                     J, edxUV, fv.number());
01081 
01082           for(MeshType::const_iterator icell(this->__mesh);
01083               not(icell.end()); ++icell) {
01084             const CellType& K = *icell;
01085 
01086             ConformTransformation T(K);
01087 
01088             TinyVector<3> massCenter;
01089             T.value(FiniteElement::massCenter(), massCenter);
01090             const real_t fValue = F(massCenter);
01091 
01092             const size_t cellNumber = icell.number();
01093 
01094             TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
01095             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
01096               index[l] = dofPositions(cellNumber,l);
01097             }
01098 
01099             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
01100               const size_t dof
01101                 = this->__degreeOfFreedomSet(testFunctionNumber,index[l]);
01102 
01103               for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
01104                 b[dof] += edxUV(l,m)*gValues[index[m]]*fValue;
01105               }
01106             }
01107           }
01108           break;
01109         }
01110         case VariationalLinearOperator::FdxV: {
01111           const VariationalOperatorFdxV& fv
01112             = dynamic_cast<const VariationalOperatorFdxV&>(*(*i));
01113 
01114           const ScalarFunctionBase& F = fv.f();
01115           const size_t testFunctionNumber = fv.testFunctionNumber();
01116 
01117           Vector<real_t> fValues (dofPositions.number());
01118 
01119           for (size_t i=0; i<dofPositions.number(); ++i) {
01120             fValues[i] = F(dofPositions.vertex(i));
01121           }
01122 
01123           MeshType::const_iterator icell0(this->__mesh);
01124           ConformTransformation T0(*icell0);
01125           JacobianTransformation J(T0);
01126 
01127           ElementaryMatrixType eUdxV = 0;
01128           generatesElementaryMatrix(PDEOperator::firstorderopTransposed,
01129                                     J, eUdxV, fv.number());
01130 
01131           for(MeshType::const_iterator icell(this->__mesh);
01132               not(icell.end()); ++icell) {
01133             const CellType& K = *icell;
01134 
01135             ConformTransformation T(K);
01136 
01137             const size_t cellNumber = icell.number();
01138 
01139             TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
01140             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
01141               index[l] = dofPositions(cellNumber,l);
01142             }
01143 
01144             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
01145               const size_t dof
01146                 = this->__degreeOfFreedomSet(testFunctionNumber,index[l]);
01147 
01148               for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
01149                 b[dof] += eUdxV(l,m)*fValues[index[m]];
01150               }
01151             }
01152           }
01153           break;
01154         }
01155         case VariationalLinearOperator::FgradGgradV: {
01156           const VariationalOperatorFgradGgradV& fv
01157             = dynamic_cast<const VariationalOperatorFgradGgradV&>(*(*i));
01158 
01159           const ScalarFunctionBase& F = fv.f();
01160           const ScalarFunctionBase& G = fv.g();
01161 
01162           const size_t testFunctionNumber = fv.testFunctionNumber();
01163 
01164           if (G.type() == ScalarFunctionBase::femfunction) {
01165             const FEMFunctionBase& gFEM = dynamic_cast<const FEMFunctionBase&>(G);
01166             if (gFEM.baseMesh() != & this->__mesh) {
01167               if (gFEM.baseMesh()->type() == Mesh::cartesianHexahedraMesh) {
01168                 const MeshType& quadratureMesh
01169                   = dynamic_cast<const MeshType&>(*gFEM.baseMesh());
01170                 for (MeshType::const_iterator iQuadratureCell(quadratureMesh);
01171                      not(iQuadratureCell.end()); ++iQuadratureCell) {
01172                   const CellType& KQuadrature = *iQuadratureCell;
01173                   ConformTransformation TQuadrature(KQuadrature);
01174                   
01175                   for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
01176                     TinyVector<3, real_t> quadrature;
01177                     TQuadrature.value(referenceQuadrature[l], quadrature);
01178 
01179                     MeshType::const_iterator icell = this->__mesh.find(quadrature);
01180                     if (icell.end()) continue;
01181 
01182                     const CellType& K = *icell;
01183                     if (not K.isFictitious()) {
01184                       const TinyVector<3,real_t> gradGValue = gFEM.gradient(quadrature);
01185                       const real_t fValue = F(quadrature);
01186 
01187                       ConformTransformation T(K);
01188                       JacobianTransformation J(T);
01189 
01190                       TinyVector<3,real_t> Xhat;
01191                       T.invertT(quadrature,Xhat);
01192 
01193                       const size_t cellNumber = icell.number();
01194 
01195                       TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
01196                       for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
01197                         index[m] = dofPositions(cellNumber,m);
01198                       }
01199 
01200                       for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
01201                         TinyVector<3,real_t> gradV;
01202 
01203                         gradV[0] = FiniteElement::instance().dxW(m,Xhat)
01204                           * J.invJacobian(0,0);
01205                         gradV[1] = FiniteElement::instance().dyW(m,Xhat)
01206                           * J.invJacobian(1,1);
01207                         gradV[2] = FiniteElement::instance().dzW(m,Xhat)
01208                           * J.invJacobian(2,2);
01209 
01210                         const size_t dof
01211                           = this->__degreeOfFreedomSet(testFunctionNumber,index[m]);
01212 
01213                         b[dof] += KQuadrature.volume()*QuadratureType::instance().weight(l)*fValue*(gradGValue*gradV);
01214                       }
01215                     }
01216                   }
01217                 }
01218                 break;
01219               }
01220             }
01221           }
01222 
01223           Vector<real_t> gValues (dofPositions.number());
01224 
01225           for (size_t i=0; i<dofPositions.number(); ++i) {
01226             gValues[i] = G(dofPositions.vertex(i));
01227           }
01228 
01229           MeshType::const_iterator icell0(this->__mesh);
01230           ConformTransformation T0(*icell0);
01231           JacobianTransformation J(T0);
01232 
01233           ElementaryMatrixType e = 0;
01234           generatesElementaryMatrix(PDEOperator::divmugrad,
01235                                     J, e);
01236 
01237           for(MeshType::const_iterator icell(this->__mesh);
01238               not(icell.end()); ++icell) {
01239             const CellType& K = *icell;
01240             ConformTransformation T(K);
01241 
01242             TinyVector<3> massCenter;
01243             T.value(FiniteElement::massCenter(), massCenter);
01244 
01245             const real_t fValue = F(massCenter);
01246 
01247             const size_t cellNumber = icell.number();
01248 
01249             TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
01250             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
01251               index[l] = dofPositions(cellNumber,l);
01252             }
01253 
01254             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
01255               const size_t dof
01256                 = __degreeOfFreedomSet(testFunctionNumber,index[l]);
01257 
01258               for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
01259                 b[dof] += e(l,m)*gValues[index[m]]*fValue;
01260               }
01261             }
01262           }
01263           break;
01264         }
01265         default: {
01266           throw ErrorHandler(__FILE__,__LINE__,
01267                              "unknown variational operator type",
01268                              ErrorHandler::unexpected);
01269         }
01270         }
01271       }
01272       break;
01273     }
01274     default: {
01275       throw ErrorHandler(__FILE__,__LINE__,
01276                          "unknown problem type",
01277                          ErrorHandler::unexpected);
01278     }
01279     }
01280     ffout(2) << "- assembling second member: done";
01281     t.stop();
01282     ffout(3) << " [cost: " << t << ']';
01283     ffout(2) << '\n';
01284   }

Here is the call graph for this function:

void BaseFEMDiscretization< Structured3DMesh , TypeOfDiscretization >::generatesElementaryVector ( ElementaryVectorType &  eVector,
const JacobianTransformation &  J,
const TinyVector< FiniteElement::numberOfQuadraturePoints, real_t > &  f 
) const [inline, protected, inherited]

Generates elementary vector

Parameters:
eVector the generated elementary vector
J the jacobian of the transformation
f the function to discretize

Definition at line 114 of file BaseFEMDiscretization.hpp.

00117   {
00118     FiniteElement::instance().integrateWj(eVector,J,f);
00119   }

void BaseFEMDiscretization< Structured3DMesh , TypeOfDiscretization >::generatesElementaryMatrix ( ElementaryMatrixSet< ElementaryMatrixType > &  eSet,
const JacobianTransformation &  J 
) const [inline, protected, inherited]

Generates elementary matrices set for a given element

Parameters:
eSet the set of elementary matrices
J the jacobian of the transformation

Definition at line 128 of file BaseFEMDiscretization.hpp.

References ElementaryMatrixSet< ElementaryMatrixType >::divMuGrad(), PDEOperator::divmugrad, PDEOperator::firstorderop, ElementaryMatrixSet< ElementaryMatrixType >::firstOrderOperatorDxUV(), ElementaryMatrixSet< ElementaryMatrixType >::firstOrderOperatorUdxV(), PDEOperator::firstorderopTransposed, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::generatesElementaryMatrix(), ElementaryMatrixSet< ElementaryMatrixType >::isDivMuGrad(), ElementaryMatrixSet< ElementaryMatrixType >::isFirstOrderDxUV(), ElementaryMatrixSet< ElementaryMatrixType >::isFirstOrderOperator(), ElementaryMatrixSet< ElementaryMatrixType >::isFirstOrderUdxV(), ElementaryMatrixSet< ElementaryMatrixType >::isMassOperator(), ElementaryMatrixSet< ElementaryMatrixType >::isSecondOrderOperator(), PDEOperator::massop, ElementaryMatrixSet< ElementaryMatrixType >::massOperator(), PDEOperator::secondorderop, and ElementaryMatrixSet< ElementaryMatrixType >::secondOrderOperator().

00130   {
00131     if (eSet.isMassOperator()) {
00132       generatesElementaryMatrix(PDEOperator::massop,
00133                                 J, eSet.massOperator());
00134     }
00135     if (eSet.isFirstOrderOperator()) {
00136       for (size_t i=0; i<3; ++i) {
00137         if (eSet.isFirstOrderUdxV(i)) {
00138           generatesElementaryMatrix(PDEOperator::firstorderopTransposed,
00139                                     J,eSet.firstOrderOperatorUdxV(i),i);
00140         }
00141         if (eSet.isFirstOrderDxUV(i)) {
00142           generatesElementaryMatrix(PDEOperator::firstorderop, J,
00143                                     eSet.firstOrderOperatorDxUV(i),i);
00144         }
00145       }
00146     }
00147     if (eSet.isSecondOrderOperator()) {
00148       for (size_t i=0; i<3; ++i)
00149         for (size_t j=0; j<3; ++j) {
00150           if (eSet.isSecondOrderOperator(i,j)) {
00151             generatesElementaryMatrix(PDEOperator::secondorderop, J,
00152                                       eSet.secondOrderOperator(i,j),i,j);
00153           }
00154         }
00155     }
00156     if (eSet.isDivMuGrad()) {
00157       generatesElementaryMatrix(PDEOperator::divmugrad, J, eSet.divMuGrad());
00158     }
00159   }

Here is the call graph for this function:

void BaseFEMDiscretization< Structured3DMesh , TypeOfDiscretization >::generatesElementaryMatrix ( const PDEOperator::Type  operatorType,
const JacobianTransformation &  J,
ElementaryMatrixType &  matelem,
const size_t  i = 0,
const size_t  j = 0 
) const [inline, protected, inherited]

Generates an elementary matrix for a given operator in an element. The row and column number can be specified when operator is not scalar: $ \partial_{x_i}(w_l)\partial_{x_j}(w_k)$ for instance.

Parameters:
operatorType type of the operator
J jacobian of the transformation
matelem generated elementary matrix
i row number
j column number

Definition at line 173 of file BaseFEMDiscretization.hpp.

References PDEOperator::divmugrad, PDEOperator::firstorderop, PDEOperator::firstorderopTransposed, PDEOperator::massop, PDEOperator::secondorderop, and ErrorHandler::unexpected.

00177   {
00178     matelem = 0;
00179     switch(operatorType) {
00180     case PDEOperator::firstorderop: {
00181       FiniteElement::instance().integrateDWjWi(matelem,i,J);
00182       matelem *= J.jacobianDet();
00183       break;
00184     }
00185     case PDEOperator::firstorderopTransposed: {
00186       FiniteElement::instance().integrateWjDWi(matelem,i,J);
00187       matelem *= J.jacobianDet();
00188       break;
00189     }
00190     case PDEOperator::divmugrad: {
00191       FiniteElement::instance().integrateDWjDWi(matelem,0,0,J);
00192       FiniteElement::instance().integrateDWjDWi(matelem,1,1,J);
00193       FiniteElement::instance().integrateDWjDWi(matelem,2,2,J);
00194       matelem *= J.jacobianDet();
00195       break;
00196     }
00197     case PDEOperator::secondorderop: {
00198       FiniteElement::instance().integrateDWjDWi(matelem,i,j,J);
00199       matelem *= J.jacobianDet();
00200       break;
00201     }
00202     case PDEOperator::massop: {
00203       FiniteElement::instance().integrateWjWi(matelem,J);
00204       matelem *= J.jacobianDet();
00205       break;
00206     }
00207     default: {
00208       throw ErrorHandler(__FILE__,__LINE__,
00209                          "unable to built elementary matrix for this operator",
00210                          ErrorHandler::unexpected);
00211     }
00212     }
00213   }

void Discretization::setDirichletList ( const Vector< bool > &  dirichletList  )  [inline, inherited]

Sets dirichlet vertices list

Parameters:
dirichletList list of dirichlet vertices

Definition at line 63 of file Discretization.hpp.

References Discretization::__dirichletList, and ASSERT.

00064   {
00065     ASSERT(__dirichletList == 0);
00066     __dirichletList = &dirichletList;
00067   }

const Problem& Discretization::problem (  )  const [inline, inherited]

Read only access to the problem

Returns:
__problem

Definition at line 109 of file Discretization.hpp.

References Discretization::__problem.

Referenced by SpectralLegendreDiscretizationNonConform::_getLocalDiscretizer(), assembleSecondMember(), and FEMDiscretization< GivenMeshType, TypeOfDiscretization >::assembleSecondMember().

00110   {
00111     return __problem;
00112   }

BaseMatrix& Discretization::A (  )  [inline, inherited]

Read only access to the discretization type

Returns:
__type Access to the matrix

__A

Definition at line 130 of file Discretization.hpp.

References Discretization::__A.

Referenced by SpectralLegendreDiscretizationNonConform::assembleMatrix(), SpectralLegendreDiscretizationConform::assembleMatrix(), assembleMatrix(), and FEMDiscretization< GivenMeshType, TypeOfDiscretization >::assembleMatrix().

00131   {
00132     return __A;
00133   }

BaseVector& Discretization::b (  )  [inline, inherited]

Access to the second member

Returns:
__b

Definition at line 140 of file Discretization.hpp.

References Discretization::__b.

Referenced by SpectralLegendreDiscretizationNonConform::_getLocalDiscretizer(), SpectralLegendreDiscretizationNonConform::assembleSecondMember(), SpectralLegendreDiscretizationConform::assembleSecondMember(), assembleSecondMember(), and FEMDiscretization< GivenMeshType, TypeOfDiscretization >::assembleSecondMember().

00141   {
00142     return __b;
00143   }

Member Data Documentation

const MeshType& BaseFEMDiscretization< Structured3DMesh , TypeOfDiscretization >::__mesh [protected, inherited]

Mesh used to perform discretization.

Definition at line 95 of file BaseFEMDiscretization.hpp.

ElementaryMatrixSet<ElementaryMatrixType> BaseFEMDiscretization< Structured3DMesh , TypeOfDiscretization >::__eSet [mutable, protected, inherited]

Set of elementary matrices.

Definition at line 98 of file BaseFEMDiscretization.hpp.

DiscretizedOperators<ElementaryMatrixType> BaseFEMDiscretization< Structured3DMesh , TypeOfDiscretization >::__discretizedOperators [mutable, protected, inherited]

Operators that are discretized.

Definition at line 101 of file BaseFEMDiscretization.hpp.

const DegreeOfFreedomSet& BaseFEMDiscretization< Structured3DMesh , TypeOfDiscretization >::__degreeOfFreedomSet [protected, inherited]

Set of degrees of freedom.

Definition at line 104 of file BaseFEMDiscretization.hpp.

const DiscretizationType Discretization::__discretizationType [protected, inherited]

The full description of the discretization

Reimplemented in SpectralLegendreDiscretizationNonConform.

Definition at line 43 of file Discretization.hpp.

Referenced by FEMDiscretization(), and FEMDiscretization< GivenMeshType, TypeOfDiscretization >::FEMDiscretization().

const Problem& Discretization::__problem [protected, inherited]

The PDEProblem to discretize.

Definition at line 46 of file Discretization.hpp.

Referenced by Discretization::problem().

BaseMatrix& Discretization::__A [protected, inherited]

The matrix which will contain the discretization.

Definition at line 49 of file Discretization.hpp.

Referenced by SpectralLegendreDiscretizationNonConform::_getLocalDiscretizer(), Discretization::A(), SpectralLegendreDiscretizationNonConform::assembleMatrix(), SpectralLegendreDiscretizationConform::assembleMatrix(), assembleMatrix(), and FEMDiscretization< GivenMeshType, TypeOfDiscretization >::assembleMatrix().

BaseVector& Discretization::__b [protected, inherited]

The second member.

Definition at line 52 of file Discretization.hpp.

Referenced by SpectralLegendreDiscretizationNonConform::_getLocalDiscretizer(), SpectralLegendreDiscretizationNonConform::assembleSecondMember(), SpectralLegendreDiscretizationConform::assembleSecondMember(), assembleSecondMember(), FEMDiscretization< GivenMeshType, TypeOfDiscretization >::assembleSecondMember(), and Discretization::b().

const Vector<bool>* Discretization::__dirichletList [protected, inherited]

elimination dirichlet informations

Definition at line 55 of file Discretization.hpp.

Referenced by Discretization::setDirichletList(), timesX(), FEMDiscretization< GivenMeshType, TypeOfDiscretization >::timesX(), transposedTimesX(), and FEMDiscretization< GivenMeshType, TypeOfDiscretization >::transposedTimesX().

The documentation for this class was generated from the following file:
Generated on Wed Nov 19 00:06:23 2008 for FreeFEM3D (aka ff3d) by  doxygen 1.5.6