FEMDiscretization< GivenMeshType, TypeOfDiscretization > Class Template Reference

#include <FEMDiscretization.hpp>

Inheritance diagram for FEMDiscretization< GivenMeshType, TypeOfDiscretization >:

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

Collaboration graph
[legend]

List of all members.

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 MeshType &m, const DiscretizationType &discretizationType, BaseMatrix &a, BaseVector &bb, const DegreeOfFreedomSet &dof)
 ~FEMDiscretization ()
void setDirichletList (const Vector< bool > &dirichletList)
const Problemproblem () const
BaseMatrixA ()
BaseVectorb ()

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

Private Types

typedef GivenMeshType 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.


Detailed Description

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
class FEMDiscretization< GivenMeshType, TypeOfDiscretization >

Definition at line 59 of file FEMDiscretization.hpp.


Member Typedef Documentation

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
typedef GivenMeshType FEMDiscretization< GivenMeshType, TypeOfDiscretization >::MeshType [private]

The type of mesh used for discretization.

Reimplemented from BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >.

Definition at line 65 of file FEMDiscretization.hpp.

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
typedef MeshType::CellType FEMDiscretization< GivenMeshType, TypeOfDiscretization >::CellType [private]

The geometry of finite elements.

Reimplemented from BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >.

Definition at line 68 of file FEMDiscretization.hpp.

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
typedef FiniteElementTraits<CellType, TypeOfDiscretization>::Type FEMDiscretization< GivenMeshType, TypeOfDiscretization >::FiniteElement [private]

The finite element.

Reimplemented from BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >.

Definition at line 72 of file FEMDiscretization.hpp.

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
typedef FiniteElement::QuadratureType FEMDiscretization< GivenMeshType, TypeOfDiscretization >::QuadratureType [private]

quadrature type

Definition at line 75 of file FEMDiscretization.hpp.

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
typedef FiniteElement::ElementaryMatrix FEMDiscretization< GivenMeshType, TypeOfDiscretization >::ElementaryMatrixType [private]

Elementary matrices type.

Reimplemented from BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >.

Definition at line 78 of file FEMDiscretization.hpp.

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
typedef FiniteElement::ElementaryVector FEMDiscretization< GivenMeshType, TypeOfDiscretization >::ElementaryVectorType [private]

Elementary vector type.

Reimplemented from BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >.

Definition at line 81 of file FEMDiscretization.hpp.

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
typedef FiniteElementTraits<CellType, TypeOfDiscretization>::Transformation FEMDiscretization< GivenMeshType, TypeOfDiscretization >::ConformTransformation [private]

type of transformation from reference element

Reimplemented from BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >.

Definition at line 85 of file FEMDiscretization.hpp.

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
typedef FiniteElementTraits<CellType, TypeOfDiscretization>::JacobianTransformation FEMDiscretization< GivenMeshType, TypeOfDiscretization >::JacobianTransformation [private]

Associated jacobian.

Reimplemented from BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >.

Definition at line 89 of file FEMDiscretization.hpp.


Constructor & Destructor Documentation

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
FEMDiscretization< GivenMeshType, TypeOfDiscretization >::FEMDiscretization ( const Problem p,
const MeshType m,
const DiscretizationType discretizationType,
BaseMatrix a,
BaseVector bb,
const DegreeOfFreedomSet dof 
) [inline]

Constructor of the discretization

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

Definition at line 623 of file FEMDiscretization.hpp.

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

00629     : BaseFEMDiscretization<MeshType,
00630                             TypeOfDiscretization>(p, m, discretizationType, a, bb, dof)
00631   {
00632     SolverInformationCenter::instance().pushMesh(&m);
00633     SolverInformationCenter::instance().pushDiscretizationType(&(this->__discretizationType));
00634   }

Here is the call graph for this function:

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
FEMDiscretization< GivenMeshType, TypeOfDiscretization >::~FEMDiscretization (  )  [inline]

Virtual destructor

Definition at line 640 of file FEMDiscretization.hpp.

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

00641   {
00642     SolverInformationCenter::instance().pop();
00643   }

Here is the call graph for this function:


Member Function Documentation

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
void FEMDiscretization< GivenMeshType, TypeOfDiscretization >::assembleMatrix (  )  [inline, virtual]

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

Implements Discretization.

Definition at line 97 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::unAssembled, and ErrorHandler::unexpected.

00098   {
00099     const ScalarDegreeOfFreedomPositionsSet& dofPositions = this->__degreeOfFreedomSet.positionsSet(0);
00100   
00101     switch ((this->__A).type()) {
00102     case BaseMatrix::doubleHashedMatrix: {
00103       Timer t;
00104       t.start();
00105 
00106       ffout(2) << "- assembling matrix\n";
00107 
00108       DoubleHashedMatrix& A = dynamic_cast<DoubleHashedMatrix&>(this->__A);
00109 
00110       for(typename MeshType::const_iterator icell(this->__mesh);
00111           not(icell.end()); ++icell) {
00112         const CellType& K = *icell;
00113 
00114         const size_t cellNumber = icell.number();
00115         TinyVector<FiniteElement::numberOfDegreesOfFreedom,
00116                    size_t> index;
00117 
00118         for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00119           index[l] = dofPositions(cellNumber,l);
00120         }
00121 
00122         ConformTransformation T(K);
00123         JacobianTransformation J(T);
00124         generatesElementaryMatrix(this->__eSet,J);
00125 
00126         TinyVector<3> massCenter;
00127         T.value(FiniteElement::massCenter(), massCenter);
00128 
00129         for(typename DiscretizedOperators<ElementaryMatrixType>
00130               ::iterator iOperator = this->__discretizedOperators.begin();
00131             iOperator != this->__discretizedOperators.end(); ++iOperator) {
00132           const real_t coef = ((*iOperator).second)(massCenter);
00133 
00134           const ElementaryMatrixType& elementaryMatrix
00135             = (*(*iOperator).first);
00136 
00137           for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
00138             const size_t I // Index of Ith value
00139               = __degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
00140             for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
00141               const size_t J // Index of Jth value
00142                 = __degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
00143               A(I,J) += coef*elementaryMatrix(l,m);
00144             }
00145           }
00146         }
00147       }
00148 
00149       ffout(2) << "- assembling matrix: done";
00150       t.stop();
00151       ffout(3) << " [cost: " << t << ']';
00152       ffout(2) << '\n';
00153       break;
00154     }
00155     case BaseMatrix::unAssembled: {
00156       UnAssembledMatrix& A = dynamic_cast<UnAssembledMatrix&>(this->__A);
00157       A.setDiscretization(this);
00158       break;
00159     }
00160 
00161     default: {
00162       throw ErrorHandler(__FILE__,__LINE__,
00163                          "unexpected matrix type",
00164                        ErrorHandler::unexpected);
00165     }
00166     }
00167   }

Here is the call graph for this function:

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
void FEMDiscretization< GivenMeshType, 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 175 of file FEMDiscretization.hpp.

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

00176   {
00177     const Vector<real_t>& x = dynamic_cast<const Vector<real_t>&>(X);
00178     Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
00179     z=0;
00180 
00181     const ScalarDegreeOfFreedomPositionsSet& dofPositions
00182       = this->__degreeOfFreedomSet.positionsSet(0);
00183   
00184     for(typename MeshType::const_iterator icell(this->__mesh);
00185         not(icell.end()); ++icell) {
00186       const CellType& K = *icell;
00187 
00188       ConformTransformation T(K);
00189       JacobianTransformation J(T);
00190       generatesElementaryMatrix(this->__eSet,J);
00191 
00192       const size_t cellNumber = icell.number();
00193 
00194       TinyVector<FiniteElement::numberOfDegreesOfFreedom,size_t> index;
00195       for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00196         index[l] = dofPositions(cellNumber,l);
00197       }
00198 
00199       TinyVector<3> massCenter;
00200       T.value(FiniteElement::massCenter(), massCenter);
00201 
00202       for(typename DiscretizedOperators<ElementaryMatrixType>
00203             ::iterator iOperator = this->__discretizedOperators.begin();
00204           iOperator != this->__discretizedOperators.end(); ++iOperator) {
00205         const real_t coef = ((*iOperator).second)(massCenter);
00206         const ElementaryMatrixType& elementaryMatrix
00207           = (*(*iOperator).first);
00208 
00209         for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
00210           const size_t I1 // Index of Ith value
00211             = __degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
00212           if (not((*this->__dirichletList)[I1])) {
00213             for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
00214               const size_t I2 // Index of Jth value
00215                 = __degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
00216               if (not((*this->__dirichletList)[I2])) {
00217                 z[I1] += coef*elementaryMatrix(l,m)*x[I2];
00218               }
00219             }
00220           }
00221         }
00222       }
00223     }
00224   }

Here is the call graph for this function:

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
void FEMDiscretization< GivenMeshType, 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 232 of file FEMDiscretization.hpp.

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

00233   {
00234     const Vector<real_t>& x = dynamic_cast<const Vector<real_t>&>(X);
00235     Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
00236     z=0;
00237 
00238     const ScalarDegreeOfFreedomPositionsSet& dofPositions
00239       = this->__degreeOfFreedomSet.positionsSet(0);
00240   
00241     for(typename MeshType::const_iterator icell(this->__mesh);
00242         not(icell.end()); ++icell) {
00243       const CellType& K = *icell;
00244 
00245       ConformTransformation T(K);
00246       JacobianTransformation J(T);
00247       generatesElementaryMatrix(this->__eSet,J);
00248 
00249       const size_t cellNumber = icell.number();
00250 
00251       TinyVector<FiniteElement::numberOfDegreesOfFreedom,size_t> index;
00252       for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00253         index[l] = dofPositions(cellNumber,l);
00254       }
00255 
00256       TinyVector<3> massCenter;
00257       T.value(FiniteElement::massCenter(), massCenter);
00258 
00259       for(typename DiscretizedOperators<ElementaryMatrixType>
00260             ::iterator iOperator = this->__discretizedOperators.begin();
00261           iOperator != this->__discretizedOperators.end(); ++iOperator) {
00262         const real_t coef = ((*iOperator).second)(massCenter);
00263         const ElementaryMatrixType& elementaryMatrix
00264           = (*(*iOperator).first);
00265 
00266         for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
00267           const size_t I1 // Index of Ith value
00268             = __degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
00269           if (not((*this->__dirichletList)[I1])) {
00270             for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
00271               const size_t I2 // Index of Jth value
00272                 = __degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
00273               if (not((*this->__dirichletList)[I2])) {
00274                 z[I2] += coef*elementaryMatrix(l,m)*x[I1];
00275               }
00276             }
00277           }
00278         }
00279       }
00280     }
00281   }

Here is the call graph for this function:

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
void FEMDiscretization< GivenMeshType, TypeOfDiscretization >::getDiagonal ( BaseVector Z  )  const [inline, virtual]

Computes diagonal of the operator

Parameters:
Z diagonal of the operator

Implements Discretization.

Definition at line 288 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().

00289   {
00290     Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
00291 
00292     const ScalarDegreeOfFreedomPositionsSet& dofPositions
00293       = this->__degreeOfFreedomSet.positionsSet(0);
00294   
00295     for(typename MeshType::const_iterator icell(this->__mesh);
00296         not(icell.end()); ++icell) {
00297       const CellType& K = *icell;
00298       const size_t cellNumber = icell.number();
00299 
00300       ConformTransformation T(K);
00301       JacobianTransformation J(T);
00302       generatesElementaryMatrix(this->__eSet,J);
00303 
00304       size_t index[FiniteElement::numberOfDegreesOfFreedom];
00305       for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00306         index[l] = dofPositions(cellNumber,l);
00307       }
00308 
00309       TinyVector<3> massCenter;
00310       T.value(FiniteElement::massCenter(), massCenter);
00311 
00312       for(typename DiscretizedOperators<ElementaryMatrixType>
00313             ::iterator iOperator = this->__discretizedOperators.begin();
00314           iOperator != this->__discretizedOperators.end(); ++iOperator) {
00315         // Only diagonal opertors have to be considerate
00316         if (((*iOperator).second).i() == ((*iOperator).second).j()) {
00317           const real_t coef = ((*iOperator).second)(massCenter);
00318 
00319           const ElementaryMatrixType& elementaryMatrix
00320             = (*(*iOperator).first);
00321 
00322           for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
00323             const size_t I // Index of Ith value
00324               = this->__degreeOfFreedomSet(((*iOperator).second).i(),
00325                                            index[l]);
00326             z[I] += coef*elementaryMatrix(l,l);
00327           }
00328         }
00329       }
00330     }
00331   }

Here is the call graph for this function:

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
void FEMDiscretization< GivenMeshType, TypeOfDiscretization >::assembleSecondMember (  )  [inline, virtual]

Second member assembling

The elementary vector

Implements Discretization.

Definition at line 337 of file FEMDiscretization.hpp.

References Discretization::__b, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__degreeOfFreedomSet, BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__mesh, Discretization::b(), VariationalProblem::beginLinearOperator(), PDEOperator::divmugrad, VariationalProblem::endLinearOperator(), VariationalOperatorFgradGgradV::f(), VariationalOperatorFdxV::f(), VariationalOperatorFdxGV::f(), VariationalOperatorFV::f(), VariationalLinearOperator::FdxGV, VariationalLinearOperator::FdxV, ffout(), VariationalLinearOperator::FgradGgradV, PDEOperator::firstorderop, PDEOperator::firstorderopTransposed, VariationalLinearOperator::FV, VariationalOperatorFgradGgradV::g(), VariationalOperatorFdxGV::g(), BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::generatesElementaryMatrix(), BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::generatesElementaryVector(), VariationalProblem::hasJumpOrMean(), ErrorHandler::normal, VariationalOperatorFdxV::number(), VariationalOperatorFdxGV::number(), ScalarDegreeOfFreedomPositionsSet::number(), Problem::numberOfUnknown(), Problem::pde, DegreeOfFreedomSet::positionsSet(), Discretization::problem(), PDESystem::secondMember(), Timer::start(), Timer::stop(), VariationalOperator::testFunctionNumber(), ErrorHandler::unexpected, Problem::variational, and ScalarDegreeOfFreedomPositionsSet::vertex().

00338   {
00339     Timer t;
00340     t.start();
00341 
00342     ffout(2) << "- assembling second member\n";
00344     ElementaryVectorType eVect;
00345     Vector<real_t>& b = (static_cast<Vector<real_t>&>(this->__b));
00346     b = 0;
00347 
00348     const QuadratureType& referenceQuadrature = QuadratureType::instance();
00349 
00350     const ScalarDegreeOfFreedomPositionsSet& dofPositions
00351       = this->__degreeOfFreedomSet.positionsSet(0);
00352   
00353     switch(this->problem().type()) {
00354     case Problem::pde: {
00355       const PDESystem& pdeSystem
00356         = static_cast<const PDESystem&>(this->problem());
00357 
00358       for(typename MeshType::const_iterator icell(this->__mesh);
00359           not(icell.end()); ++icell) {
00360         const CellType& K = *icell;
00361         eVect = 0;
00362         ConformTransformation T(K);
00363         JacobianTransformation J(T);
00364 
00365         TinyVector<QuadratureType::numberOfQuadraturePoints,
00366                    TinyVector<3> > quadrature;
00367         for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
00368           T.value(referenceQuadrature[l], quadrature[l]);
00369         }
00370 
00371         const size_t cellNumber = icell.number();
00372 
00373         TinyVector<FiniteElement::numberOfDegreesOfFreedom, int> index;
00374         for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00375           index[l] = dofPositions(cellNumber,l);
00376         }
00377 
00378         for (size_t i=0; i<this->problem().numberOfUnknown(); ++i) {
00379           TinyVector<FiniteElement::numberOfQuadraturePoints, real_t> f;
00380           const ScalarFunctionBase& F = pdeSystem.secondMember(i);
00381 
00382           for (size_t l=0; l<FiniteElement::numberOfQuadraturePoints; ++l)
00383             f[l] = F(quadrature[l]);
00384 
00385           generatesElementaryVector(eVect,J,f);
00386 
00387           for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00388             const size_t dof = __degreeOfFreedomSet(i,index[l]);
00389             b[dof] += eVect[l]*J.jacobianDet();
00390           }
00391         }
00392       }
00393       break;
00394     }
00395     case Problem::variational: {
00396       const VariationalProblem& variationalProblem
00397         =  dynamic_cast<const VariationalProblem&>(this->problem());
00398 
00399       if (variationalProblem.hasJumpOrMean()) {
00400         throw ErrorHandler(__FILE__,__LINE__,
00401                            "cannot discretize problems with jumps or means in FEM !",
00402                            ErrorHandler::normal);
00403       }
00404 
00405       for (VariationalProblem::linearOperatorConst_iterator i
00406              = variationalProblem.beginLinearOperator();
00407            i != variationalProblem.endLinearOperator(); ++i) {
00408         switch ((*(*i)).type()) {
00409         case VariationalLinearOperator::FV: {
00410           const VariationalOperatorFV& fv
00411             = dynamic_cast<const VariationalOperatorFV&>(*(*i));
00412 
00413           const ScalarFunctionBase& F = fv.f();
00414           const size_t testFunctionNumber = fv.testFunctionNumber();
00415 
00416           for(typename MeshType::const_iterator icell(this->__mesh);
00417               not(icell.end()); ++icell) {
00418             const CellType& K = *icell;
00419             eVect = 0;
00420             ConformTransformation T(K);
00421             JacobianTransformation J(T);
00422 
00423             TinyVector<QuadratureType::numberOfQuadraturePoints,
00424                        TinyVector<3> > quadrature;
00425             for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
00426               T.value(referenceQuadrature[l], quadrature[l]);
00427             }
00428 
00429             const size_t cellNumber = icell.number();
00430 
00431             TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
00432             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00433               index[l] = dofPositions(cellNumber,l);
00434             }
00435 
00436             TinyVector<FiniteElement::numberOfQuadraturePoints,real_t> f;
00437 
00438             for (size_t l=0; l<FiniteElement::numberOfQuadraturePoints; ++l)
00439               f[l] = F(quadrature[l]);
00440 
00441             generatesElementaryVector(eVect,J,f);
00442 
00443             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00444               const size_t dof = __degreeOfFreedomSet(testFunctionNumber,index[l]);
00445               b[dof] += eVect[l]*J.jacobianDet();
00446             }
00447           }
00448           break;
00449         }
00450         case VariationalLinearOperator::FdxGV: {
00451           const VariationalOperatorFdxGV& fv
00452             = dynamic_cast<const VariationalOperatorFdxGV&>(*(*i));
00453 
00454           const ScalarFunctionBase& F = fv.f();
00455           const ScalarFunctionBase& G = fv.g();
00456           const size_t testFunctionNumber = fv.testFunctionNumber();
00457 
00458           Vector<real_t> gValues (dofPositions.number());
00459 
00460           for (size_t i=0; i<dofPositions.number(); ++i) {
00461             gValues[i] = G(dofPositions.vertex(i));
00462           }
00463 
00464 
00465           for (typename MeshType::const_iterator icell(this->__mesh);
00466                not(icell.end()); ++icell) {
00467             const CellType& K = *icell;
00468 
00469             ConformTransformation T(K);
00470             JacobianTransformation J(T);
00471 
00472             ElementaryMatrixType edxUV = 0;
00473             generatesElementaryMatrix(PDEOperator::firstorderop,
00474                                       J, edxUV, fv.number());
00475 
00476             TinyVector<3> massCenter;
00477             T.value(FiniteElement::massCenter(), massCenter);
00478 
00479             const real_t fValue = F(massCenter);
00480 
00481             const size_t cellNumber = icell.number();
00482 
00483             TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
00484             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00485               index[l] = dofPositions(cellNumber,l);
00486             }
00487 
00488             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00489               const size_t dof
00490                 = __degreeOfFreedomSet(testFunctionNumber,index[l]);
00491 
00492               for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
00493                 b[dof] += edxUV(l,m)*gValues[index[m]]*fValue;
00494               }
00495             }
00496           }
00497           break;
00498         }
00499         case VariationalLinearOperator::FdxV: {
00500           const VariationalOperatorFdxV& fv
00501             = dynamic_cast<const VariationalOperatorFdxV&>(*(*i));
00502 
00503           const ScalarFunctionBase& F = fv.f();
00504           const size_t testFunctionNumber = fv.testFunctionNumber();
00505 
00506           Vector<real_t> fValues (dofPositions.number());
00507 
00508           for (size_t i=0; i<dofPositions.number(); ++i) {
00509             fValues[i] = F(dofPositions.vertex(i));
00510           }
00511 
00512           for(typename MeshType::const_iterator icell(this->__mesh);
00513               not(icell.end()); ++icell) {
00514             const CellType& K = *icell;
00515 
00516             ConformTransformation T(K);
00517             JacobianTransformation J(T);
00518 
00519             ElementaryMatrixType eUdxV = 0;
00520             generatesElementaryMatrix(PDEOperator::firstorderopTransposed,
00521                                       J, eUdxV, fv.number());
00522 
00523             const size_t cellNumber = icell.number();
00524 
00525             TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
00526             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00527               index[l] = dofPositions(cellNumber,l);
00528             }
00529 
00530             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00531               const size_t dof
00532                 = this->__degreeOfFreedomSet(testFunctionNumber,index[l]);
00533 
00534               for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
00535                 b[dof] += eUdxV(l,m)*fValues[index[m]];
00536               }
00537             }
00538           }
00539           break;
00540         }
00541         case VariationalLinearOperator::FgradGgradV: {
00542           const VariationalOperatorFgradGgradV& fv
00543             = dynamic_cast<const VariationalOperatorFgradGgradV&>(*(*i));
00544 
00545           const ScalarFunctionBase& F = fv.f();
00546           const ScalarFunctionBase& G = fv.g();
00547           const size_t testFunctionNumber = fv.testFunctionNumber();
00548 
00549           Vector<real_t> gValues (dofPositions.number());
00550 
00551           for (size_t i=0; i<dofPositions.number(); ++i) {
00552             gValues[i] = G(dofPositions.vertex(i));
00553           }
00554 
00555           for(typename MeshType::const_iterator icell(this->__mesh);
00556               not(icell.end()); ++icell) {
00557             const CellType& K = *icell;
00558 
00559             ConformTransformation T(K);
00560             JacobianTransformation J(T);
00561 
00562             ElementaryMatrixType e = 0;
00563             generatesElementaryMatrix(PDEOperator::divmugrad,
00564                                       J, e);
00565 
00566             TinyVector<3> massCenter;
00567             T.value(FiniteElement::massCenter(), massCenter);
00568 
00569             const real_t fValue = F(massCenter);
00570 
00571             const size_t cellNumber = icell.number();
00572 
00573             TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
00574             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00575               index[l] = dofPositions(cellNumber,l);
00576             }
00577 
00578             for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
00579               const size_t dof
00580                 = __degreeOfFreedomSet(testFunctionNumber,index[l]);
00581 
00582               for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
00583                 b[dof] += e(l,m)*gValues[index[m]]*fValue;
00584               }
00585             }
00586           }
00587           break;
00588         }
00589         default: {
00590           throw ErrorHandler(__FILE__,__LINE__,
00591                              "unknown variational operator",
00592                              ErrorHandler::unexpected);
00593         }
00594         }
00595       }
00596       break;
00597     }
00598     default: {
00599       throw ErrorHandler(__FILE__,__LINE__,
00600                          "unknown problem type",
00601                          ErrorHandler::unexpected);
00602     }
00603     }
00604     ffout(2) << "- assembling second member: done";
00605     t.stop();
00606     ffout(3) << " [cost: " << t << ']';
00607     ffout(2) << '\n';
00608   }

Here is the call graph for this function:

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
void BaseFEMDiscretization< GivenMeshType, 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.

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

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

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
void BaseFEMDiscretization< GivenMeshType, 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.

Referenced by FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::assembleMatrix(), FEMDiscretization< GivenMeshType, TypeOfDiscretization >::assembleMatrix(), FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::assembleSecondMember(), FEMDiscretization< GivenMeshType, TypeOfDiscretization >::assembleSecondMember(), BaseFEMDiscretization< Structured3DMesh, TypeOfDiscretization >::generatesElementaryMatrix(), FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::getDiagonal(), FEMDiscretization< GivenMeshType, TypeOfDiscretization >::getDiagonal(), FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::timesX(), FEMDiscretization< GivenMeshType, TypeOfDiscretization >::timesX(), FEMDiscretization< Structured3DMesh, TypeOfDiscretization >::transposedTimesX(), and FEMDiscretization< GivenMeshType, TypeOfDiscretization >::transposedTimesX().

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   }

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
void BaseFEMDiscretization< GivenMeshType, 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.

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]

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

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


Member Data Documentation

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
const MeshType& BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__mesh [protected, inherited]

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
ElementaryMatrixSet<ElementaryMatrixType> BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__eSet [mutable, protected, inherited]

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
DiscretizedOperators<ElementaryMatrixType> BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__discretizedOperators [mutable, protected, inherited]

template<typename GivenMeshType, ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
const DegreeOfFreedomSet& BaseFEMDiscretization< GivenMeshType, TypeOfDiscretization >::__degreeOfFreedomSet [protected, inherited]

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]

BaseVector& Discretization::__b [protected, inherited]

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


The documentation for this class was generated from the following file:

Generated on Wed Nov 19 00:06:18 2008 for FreeFEM3D (aka ff3d) by  doxygen 1.5.6