17 |
|
|
18 |
// $Id$ |
// $Id$ |
19 |
|
|
|
|
|
20 |
#ifndef FEM_DISCRETIZATION_HPP |
#ifndef FEM_DISCRETIZATION_HPP |
21 |
#define FEM_DISCRETIZATION_HPP |
#define FEM_DISCRETIZATION_HPP |
22 |
|
|
39 |
|
|
40 |
#include <ElementaryMatrixSet.hpp> |
#include <ElementaryMatrixSet.hpp> |
41 |
|
|
42 |
#include <Discretization.hpp> |
#include <BaseFEMDiscretization.hpp> |
43 |
|
|
44 |
#include <Mesh.hpp> |
#include <Mesh.hpp> |
45 |
#include <Structured3DMesh.hpp> |
#include <Structured3DMesh.hpp> |
46 |
|
|
|
#include <DiscretizedOperators.hpp> |
|
|
#include <DegreeOfFreedomSet.hpp> |
|
|
|
|
|
#include <FiniteElementTraits.hpp> |
|
|
|
|
|
#include <PDE.hpp> |
|
|
#include <PDEProblem.hpp> |
|
|
#include <MassOperator.hpp> |
|
|
#include <FirstOrderOperator.hpp> |
|
|
#include <DivMuGrad.hpp> |
|
|
#include <SecondOrderOperator.hpp> |
|
|
|
|
47 |
#include <Timer.hpp> |
#include <Timer.hpp> |
48 |
#include <Q1FiniteElement.hpp> |
#include <Q1FiniteElement.hpp> |
49 |
|
|
50 |
#include <DoubleHashedMatrix.hpp> |
#include <DoubleHashedMatrix.hpp> |
51 |
#include <UnAssembledMatrix.hpp> |
#include <UnAssembledMatrix.hpp> |
52 |
|
|
|
#include <VariationalProblem.hpp> |
|
|
|
|
|
#include <VariationalOperatorFV.hpp> |
|
|
#include <VariationalOperatorFdxGV.hpp> |
|
|
#include <VariationalOperatorFgradGgradV.hpp> |
|
|
|
|
|
#include <ConformTransformation.hpp> |
|
|
|
|
53 |
#warning Should not use language classes here |
#warning Should not use language classes here |
54 |
#include <FunctionExpression.hpp> |
#include <FunctionExpression.hpp> |
55 |
#include <MeshExpression.hpp> |
#include <MeshExpression.hpp> |
57 |
|
|
58 |
template <typename GivenMeshType> |
template <typename GivenMeshType> |
59 |
class FEMDiscretization |
class FEMDiscretization |
60 |
: public Discretization |
: public BaseFEMDiscretization<GivenMeshType> |
61 |
{ |
{ |
62 |
public: |
private: |
63 |
/// The type of mesh used for discretization |
/// The type of mesh used for discretization |
64 |
typedef GivenMeshType MeshType; |
typedef GivenMeshType MeshType; |
65 |
|
|
84 |
/// Associated jacobian |
/// Associated jacobian |
85 |
typedef typename FiniteElement::JacobianTransformation JacobianTransformation; |
typedef typename FiniteElement::JacobianTransformation JacobianTransformation; |
86 |
|
|
|
private: |
|
|
/// Mesh used to perform discretization |
|
|
MeshType& __mesh; |
|
|
|
|
|
/// Set of elementary matrices |
|
|
mutable ElementaryMatrixSet <ElementaryMatrixType> __eSet; |
|
|
|
|
|
/// Operators that are discretized |
|
|
mutable DiscretizedOperators<ElementaryMatrixType> __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<ElementaryMatrixType>& 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); |
|
|
} |
|
|
} |
|
|
} |
|
|
|
|
87 |
public: |
public: |
88 |
/** |
/** |
89 |
* Assembles the matrix associated to the PDE operators of the PDE |
* Assembles the matrix associated to the PDE operators of the PDE |
314 |
} |
} |
315 |
} |
} |
316 |
|
|
|
/** |
|
|
* Access function to the discretization mesh |
|
|
* |
|
|
* @return the mesh |
|
|
*/ |
|
|
Mesh& mesh() |
|
|
{ |
|
|
return (__mesh); |
|
|
} |
|
|
|
|
|
/** |
|
|
* Read only access to the discretization mesh |
|
|
* |
|
|
* @return the mesh |
|
|
*/ |
|
|
const Mesh& mesh() const |
|
|
{ |
|
|
return (__mesh); |
|
|
} |
|
|
|
|
|
|
|
317 |
public: |
public: |
318 |
|
|
319 |
/** |
/** |
331 |
BaseMatrix& a, |
BaseMatrix& a, |
332 |
BaseVector& bb, |
BaseVector& bb, |
333 |
const DegreeOfFreedomSet& dof) |
const DegreeOfFreedomSet& dof) |
334 |
: Discretization(Discretization::FEM,p,a,bb), |
: BaseFEMDiscretization<MeshType>(p, m, a, bb, dof) |
|
__mesh(m), |
|
|
__eSet(problem()), |
|
|
__discretizedOperators(__eSet,problem()), |
|
|
__degreeOfFreedomSet(dof) |
|
335 |
{ |
{ |
336 |
; |
; |
337 |
} |
} |
340 |
* Virtual destructor |
* Virtual destructor |
341 |
* |
* |
342 |
*/ |
*/ |
343 |
virtual ~FEMDiscretization() |
~FEMDiscretization() |
344 |
{ |
{ |
345 |
; |
; |
346 |
} |
} |
354 |
*/ |
*/ |
355 |
template <> |
template <> |
356 |
class FEMDiscretization<Structured3DMesh> |
class FEMDiscretization<Structured3DMesh> |
357 |
: public Discretization |
: public BaseFEMDiscretization<Structured3DMesh> |
358 |
{ |
{ |
359 |
public: |
public: |
360 |
/// The type of mesh used for discretization |
/// The type of mesh used for discretization |
381 |
/// Associated jacobian |
/// Associated jacobian |
382 |
typedef FiniteElement::JacobianTransformation JacobianTransformation; |
typedef FiniteElement::JacobianTransformation JacobianTransformation; |
383 |
|
|
|
private: |
|
|
/// Mesh used to perform discretization |
|
|
MeshType& __mesh; |
|
|
|
|
|
/// Set of elementary matrices |
|
|
mutable ElementaryMatrixSet <ElementaryMatrixType> __eSet; |
|
|
|
|
|
/// Operators that are discretized |
|
|
mutable DiscretizedOperators<ElementaryMatrixType> __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<TinyMatrix<8,8> >& 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); |
|
|
} |
|
|
} |
|
|
} |
|
|
|
|
384 |
public: |
public: |
385 |
/** |
/** |
386 |
* Assembles the matrix associated to the PDE operators of the PDE |
* Assembles the matrix associated to the PDE operators of the PDE |
860 |
} |
} |
861 |
} |
} |
862 |
|
|
|
/** |
|
|
* Access function to the discretization mesh |
|
|
* |
|
|
* @return the mesh |
|
|
*/ |
|
|
Mesh& mesh() |
|
|
{ |
|
|
return (__mesh); |
|
|
} |
|
|
|
|
|
/** |
|
|
* Read only access to the discretization mesh |
|
|
* |
|
|
* @return the mesh |
|
|
*/ |
|
|
const Mesh& mesh() const |
|
|
{ |
|
|
return (__mesh); |
|
|
} |
|
|
|
|
|
|
|
863 |
public: |
public: |
864 |
|
|
865 |
/** |
/** |
877 |
BaseMatrix& a, |
BaseMatrix& a, |
878 |
BaseVector& bb, |
BaseVector& bb, |
879 |
const DegreeOfFreedomSet& dof) |
const DegreeOfFreedomSet& dof) |
880 |
: Discretization(Discretization::FEM,p,a,bb), |
: BaseFEMDiscretization<Structured3DMesh>(p, m, a, bb, dof) |
|
__mesh(m), |
|
|
__eSet(problem()), |
|
|
__discretizedOperators(__eSet,problem()), |
|
|
__degreeOfFreedomSet(dof) |
|
881 |
{ |
{ |
882 |
; |
; |
883 |
} |
} |
884 |
|
|
885 |
/** |
/** |
886 |
* Virtual destructor |
* destructor |
887 |
* |
* |
888 |
*/ |
*/ |
889 |
virtual ~FEMDiscretization() |
~FEMDiscretization() |
890 |
{ |
{ |
891 |
; |
; |
892 |
} |
} |