// This file is part of ff3d - http://www.freefem.org/ff3d // Copyright (C) 2001, 2002, 2003 Stéphane Del Pino // This program is free software; you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation; either version 2, or (at your option) // any later version. // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // You should have received a copy of the GNU General Public License // along with this program; if not, write to the Free Software Foundation, // Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. // $Id: BaseFEMDiscretization.hpp,v 1.1 2003/04/15 17:47:28 delpinux Exp $ #ifndef BASE_FEM_DISCRETIZATION_HPP #define BASE_FEM_DISCRETIZATION_HPP #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /** * @file BaseFEMDiscretization.hpp * @author Stéphane Del Pino * @date Mon Apr 14 00:09:03 2003 * * @brief Finite element base class * * This class is used to define finite element method standard tools * [elementary matrices generation...]. It main purpose is to * factorise pieces of code. */ template class BaseFEMDiscretization : public Discretization { protected: /// The type of mesh used for discretization typedef GivenMeshType MeshType; /// The geometry of finite elements typedef typename MeshType::ElementsGeometry ElementsGeometry; /// The finite element typedef FiniteElementTraits FiniteElement; /// Finite element type typedef typename FiniteElement::Type FiniteElementType; /// Elementary matrices type typedef typename FiniteElement::ElementaryMatrix ElementaryMatrixType; /// Elementary vector type typedef typename FiniteElement::ElementaryVector ElementaryVectorType; /// type of transformation from reference element typedef typename FiniteElement::Transformation ConformTransformation; /// Associated jacobian typedef typename FiniteElement::JacobianTransformation JacobianTransformation; /// Mesh used to perform discretization MeshType& __mesh; /// Set of elementary matrices mutable ElementaryMatrixSet __eSet; /// Operators that are discretized mutable DiscretizedOperators __discretizedOperators; /// Set of degrees of freedom const DegreeOfFreedomSet& __degreeOfFreedomSet; /** * Generates elementary vector * * @param eVector the generated elementary vector * @param J the jacobian of the transformation * @param f the function to discretize */ void generatesElementaryVector(ElementaryVectorType& eVector, const JacobianTransformation& J, const ElementaryVectorType& f) const { FiniteElementType::instance().integrateWj(eVector,J,f); } /** * Generates elementary matrices set for a given element * * @param eSet the set of elementary matrices * @param J the jacobian of the transformation */ void generatesElementaryMatrix(ElementaryMatrixSet& eSet, const JacobianTransformation& J) const { if (eSet.isMassOperator()) { generatesElementaryMatrix(PDEOperator::massop, J, eSet.massOperator()); } if (eSet.isFirstOrderOperator()) { for (size_t i=0; i<3; ++i) { if (eSet.isFirstOrderUdxV(i)) { generatesElementaryMatrix(PDEOperator::firstorderopTransposed, J,eSet.firstOrderOperatorUdxV(i),i); } if (eSet.isFirstOrderDxUV(i)) { generatesElementaryMatrix(PDEOperator::firstorderop, J, eSet.firstOrderOperatorDxUV(i),i); } } } if (eSet.isSecondOrderOperator()) { for (size_t i=0; i<3; ++i) for (size_t j=0; j<3; ++j) { if (eSet.isSecondOrderOperator(i,j)) { generatesElementaryMatrix(PDEOperator::secondorderop, J, eSet.secondOrderOperator(i,j),i,j); } } } if (eSet.isDivMuGrad()) { generatesElementaryMatrix(PDEOperator::divmugrad, J, eSet.divMuGrad()); } } /** * Generates an elementary matrix for a given operator in an * element. The row and column number can be specified when operator * is not scalar: \f$ \partial x_i(w_l)\partial x_j(w_k)\f$ for instance. * * @param operatorType type of the operator * @param J jacobian of the transformation * @param matelem generated elementary matrix * @param i row number * @param j column number */ void generatesElementaryMatrix(const PDEOperator::Type operatorType, const JacobianTransformation& J, ElementaryMatrixType& matelem, const size_t i = 0, const size_t j = 0) const { matelem = 0; switch(operatorType) { case PDEOperator::firstorderop: { FiniteElementType::instance().integrateDWjWi(matelem,i,J); matelem *= J.jacobianDet(); break; } case PDEOperator::firstorderopTransposed: { FiniteElementType::instance().integrateWjDWi(matelem,i,J); matelem *= J.jacobianDet(); break; } case PDEOperator::divmugrad: { FiniteElementType::instance().integrateDWjDWi(matelem,0,0,J); FiniteElementType::instance().integrateDWjDWi(matelem,1,1,J); FiniteElementType::instance().integrateDWjDWi(matelem,2,2,J); matelem *= J.jacobianDet(); break; } case PDEOperator::secondorderop: { FiniteElementType::instance().integrateDWjDWi(matelem,i,j,J); matelem *= J.jacobianDet(); break; } case PDEOperator::massop: { FiniteElementType::instance().integrateWjWi(matelem,J); matelem *= J.jacobianDet(); break; } default: { fferr(2) << '\n' << __FILE__ << ':' << __LINE__ << ':' << "Not implemented\n"; std::exit(1); } } } /** * Constructor of the discretization * * @param p the problem * @param m the mesh used for discretization * @param a matrix storing discretization * @param bb vector that stores second member discretization * @param dof degrees of freedom set * */ BaseFEMDiscretization(const Problem& p, MeshType& m, BaseMatrix& a, BaseVector& bb, const DegreeOfFreedomSet& dof) : Discretization(Discretization::FEM, p, a, bb), __mesh(m), __eSet(problem()), __discretizedOperators(__eSet,problem()), __degreeOfFreedomSet(dof) { ; } /** * virtual destructor * */ virtual ~BaseFEMDiscretization() { ; } }; #endif // BASE_FEM_DISCRETIZATION_HPP