v0.5.86
PoissonOperators.hpp
/**
* \file PoissonOperators.hpp
* \example PoissonOperators.hpp
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM 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 Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#ifndef __POISSONOPERATORS_HPP__
#define __POISSONOPERATORS_HPP__
namespace PoissonExample {
/**
* \brief Calculate the grad-grad operator and assemble matrix
*
* Calculate
* \f[
* \mathbf{K}=\int_\Omega \nabla \boldsymbol\phi \cdot \nabla \boldsymbol\phi \textrm{d}\Omega
* \f]
* and assemble to global matrix.
*
* This operator is executed on element for each unique combination of entities.
*
*/
struct OpK: public VolumeElementForcesAndSourcesCore::UserDataOperator {
OpK(bool symm = true):
VolumeElementForcesAndSourcesCore::UserDataOperator("U","U",OPROWCOL,symm) {
}
/**
* \brief Do calculations for give operator
* @param row_side row side number (local number) of entity on element
* @param col_side column side number (local number) of entity on element
* @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
* @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
* @param row_data data for row
* @param col_data data for column
* @return error code
*/
PetscErrorCode doWork(
int row_side,int col_side,
EntityType row_type,EntityType col_type,
DataForcesAndSourcesCore::EntData &row_data,DataForcesAndSourcesCore::EntData &col_data
) {
PetscFunctionBegin;
// get number of dofs on row
nbRows = row_data.getIndices().size();
// if no dofs on row, exit that work, nothing to do here
if(!nbRows) PetscFunctionReturn(0);
// get number of dofs on column
nbCols = col_data.getIndices().size();
// if no dofs on Columbia, exit nothing to do here
if(!nbCols) PetscFunctionReturn(0);
// get number of integration points
nbIntegrationPts = getGaussPts().size2();
// chekk if entity block is on matrix diagonal
if(
row_side==col_side&&
row_type==col_type
) {
isDiag = true; // yes, it is
} else {
isDiag = false;
}
// integrate local matrix for entity block
ierr = iNtegrate(row_data,col_data); CHKERRQ(ierr);
// asseble local matrix
ierr = aSsemble(row_data,col_data); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
protected:
///< error code
int nbRows; ///< number of dofs on rows
int nbCols; ///< number if dof on column
int nbIntegrationPts; ///< number of integration points
bool isDiag; ///< true if this block is on diagonal
FTensor::Index<'i',3> i; ///< summit Index
MatrixDouble locMat; ///< local entity block matrix
/**
* \brief Integrate grad-grad operator
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
virtual PetscErrorCode iNtegrate(
DataForcesAndSourcesCore::EntData &row_data,DataForcesAndSourcesCore::EntData &col_data
) {
PetscFunctionBegin;
// set size of local entity bock
locMat.resize(nbRows,nbCols,false);
// clear matrux
locMat.clear();
// get element volume
double vol = getVolume();
// get integration weigths
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
// get base function gradient on rows
FTensor::Tensor1<double*,3> t_row_grad = row_data.getFTensor1DiffN<3>();
// loop over integration points
for(int gg = 0;gg!=nbIntegrationPts;gg++) {
// take into account Jacobean
const double alpha = t_w*vol;
// noalias(locMat) += alpha*prod(row_data.getDiffN(gg),trans(col_data.getDiffN(gg)));
// take fist element to local matrix
FTensor::Tensor0<double*> a(&*locMat.data().begin());
// loop over rows base functions
for(int rr = 0;rr!=nbRows;rr++) {
// get column base functions gradient at gauss point gg
FTensor::Tensor1<double*,3> t_col_grad = col_data.getFTensor1DiffN<3>(gg,0);
// loop over columbs
for(int cc = 0;cc!=nbCols;cc++) {
// calculate element of loacl matrix
a += alpha*(t_row_grad(i)*t_col_grad(i));
++t_col_grad; // move to another gradient of base function on column
++a; // move to another element of local matrix in column
}
++t_row_grad; // move to another element of gradient of base function on row
}
++t_w; // move to another integration weight
}
PetscFunctionReturn(0);
}
/**
* \brief Assemble local entity block matrix
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
virtual PetscErrorCode aSsemble(
DataForcesAndSourcesCore::EntData &row_data,DataForcesAndSourcesCore::EntData &col_data
) {
PetscFunctionBegin;
// get pointer to first global index on row
const int* row_indices = &*row_data.getIndices().data().begin();
// get pointer to first global index on column
const int* col_indices = &*col_data.getIndices().data().begin();
Mat B = getFEMethod()->ksp_B!=PETSC_NULL? getFEMethod()->ksp_B : getFEMethod()->snes_B;
// assemble local matrix
ierr = MatSetValues(
B, nbRows,row_indices,nbCols,col_indices,&*locMat.data().begin(),ADD_VALUES
if(!isDiag&&sYmm) {
// if not diagonal term and since global matrix is symmetric assemble
// transpose term.
locMat = trans(locMat);
ierr = MatSetValues(
B,nbCols,col_indices,nbRows,row_indices,&*locMat.data().begin(),ADD_VALUES
}
PetscFunctionReturn(0);
}
};
/**
* \brief template class for integration oh the right hand side
*/
template<typename OPBASE>
struct OpBaseRhs: public OPBASE {
OpBaseRhs(const std::string field_name):
OPBASE(field_name,OPBASE::OPROW) {
}
/**
* \brief This function is called by finite element
*
* Do work is composed from two operations, integrate and assembly. Also,
* it set values nbRows, and nbIntegrationPts.
*
*/
PetscErrorCode doWork(
int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
) {
PetscFunctionBegin;
// get number of dofs on row
nbRows = row_data.getIndices().size();
if(!nbRows) PetscFunctionReturn(0);
// get number of integration points
nbIntegrationPts = OPBASE::getGaussPts().size2();
// integrate local vector
ierr = iNtegrate(row_data); CHKERRQ(ierr);
// assemble local vector
ierr = aSsemble(row_data); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
/**
* \brief Class dedicated to integrate operator
* @param data entity data on element row
* @return error code
*/
virtual PetscErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &data) = 0;
/**
* \brief Class dedicated to assemble operator to global system vector
* @param data entity data (indices, base functions, etc. ) on element row
* @return error code
*/
virtual PetscErrorCode aSsemble(DataForcesAndSourcesCore::EntData &data) = 0;
protected:
///< error code
int nbRows; ///< number of dofs on row
int nbIntegrationPts; ///< number of integration points
};
/**
* \brief Operator calculate source term,
*
* \f[
* \mathbf{F} = \int_\Omega \boldsymbol\phi f \textrm{d}\Omega
* \f]
*
*/
struct OpF: public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
typedef boost::function<double (const double,const double,const double)> FSource;
OpF(FSource f_source):
OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator>("U"),
fSource(f_source) {
}
protected:
FSource fSource;
VectorDouble locVec;
/**
* \brief Integrate local entity vector
* @param data entity data on element row
* @return error code
*/
PetscErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &data) {
PetscFunctionBegin;
// set size of local vector
locVec.resize(nbRows,false);
// clear local entity vector
locVec.clear();
// get finite element volume
double vol = getVolume();
// get integration weights
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
// get base functions on entity
FTensor::Tensor0<double*> t_v = data.getFTensor0N();
// get coordinates at integration points
FTensor::Tensor1<double*,3> t_coords = getTensor1CoordsAtGaussPts();
// loop over all integration points
for(int gg = 0;gg!=nbIntegrationPts;gg++) {
// evaluate constant term
const double alpha = vol*t_w*fSource(t_coords(NX),t_coords(NY),t_coords(NZ));
// get element of local vector
FTensor::Tensor0<double*> t_a(&*locVec.data().begin());
// loop over base functions
for(int rr = 0;rr!=nbRows;rr++) {
// add to local vector source term
t_a -= alpha*t_v;
++t_a; // move to next element of local vector
++t_v; // move to next base function
}
++t_w; // move to next integration weight
++t_coords; // move to next physical coordinates at integration point
}
PetscFunctionReturn(0);
}
/**
* \brief assemble local entity vector to the global right hand side
* @param data entity data, i.e. global indices of local vector
* @return error code
*/
PetscErrorCode aSsemble(DataForcesAndSourcesCore::EntData &data) {
PetscFunctionBegin;
// get global indices of local vector
const int* indices = &*data.getIndices().data().begin();
// get values from local vector
const double* vals = &*locVec.data().begin();
Vec f = getFEMethod()->ksp_f!=PETSC_NULL ? getFEMethod()->ksp_f : getFEMethod()->snes_f;
// assemble vector
ierr = VecSetValues(f,nbRows,indices,vals,ADD_VALUES); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
};
/**
* \brief Calculate constrains matrix
*
* \f[
* \mathbf{C} = \int_{\partial\Omega} \boldsymbol\psi \boldsymbol\phi \textrm{d}\partial\Omega
* \f]
* where \f$\lambda \f$ is base function on boundary
*
*/
struct OpC: public FaceElementForcesAndSourcesCore::UserDataOperator {
OpC(const bool assemble_transpose):
FaceElementForcesAndSourcesCore::UserDataOperator("L","U",OPROWCOL,false),
assembleTraspose(assemble_transpose) {
}
PetscErrorCode doWork(
int row_side,int col_side,
EntityType row_type,EntityType col_type,
DataForcesAndSourcesCore::EntData &row_data,DataForcesAndSourcesCore::EntData &col_data
) {
PetscFunctionBegin;
// get number of dofs on row
nbRows = row_data.getIndices().size();
// exit here if no dofs on row, nothing to do
if(!nbRows) PetscFunctionReturn(0);
// get number of dofs on column,
nbCols = col_data.getIndices().size();
// exit here if no dofs on roe, nothing to do
if(!nbCols) PetscFunctionReturn(0);
// get number of integration points
nbIntegrationPts = getGaussPts().size2();
// integrate local constrains matrix
ierr = iNtegrate(row_data,col_data); CHKERRQ(ierr);
// assemble local constrains matrix
ierr = aSsemble(row_data,col_data); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
private:
///< error code
int nbRows; ///< number of dofs on row
int nbCols; ///< number of dofs on column
int nbIntegrationPts; ///< number of integration points
const bool assembleTraspose; ///< assemble transpose, i.e. CT if set to true
MatrixDouble locMat; ///< local constrains matrxi
/** \brief Integrate local constrains matrix
*/
inline PetscErrorCode iNtegrate(
DataForcesAndSourcesCore::EntData &row_data,DataForcesAndSourcesCore::EntData &col_data
) {
PetscFunctionBegin;
// set size of local constrains matrix
locMat.resize(nbRows,nbCols,false);
// clear matrix
locMat.clear();
// get area of element
const double area = getArea();
// get integration weights
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
// get base functions on entity
FTensor::Tensor0<double*> t_row = row_data.getFTensor0N();
// run over integration points
for(int gg = 0;gg!=nbIntegrationPts;gg++) {
const double alpha = area*t_w;
// get element of local matrix
FTensor::Tensor0<double*> c(&*locMat.data().begin());
// run over base functions on rows
for(int rr = 0;rr!=nbRows;rr++) {
// get first base functions on column for integration point gg
FTensor::Tensor0<double*> t_col = col_data.getFTensor0N(gg,0);
// run over base function on column
for(int cc = 0;cc!=nbCols;cc++) {
// integrate element of constrains matrix
c += alpha*t_row*t_col;
++t_col; // move to next base function on column
++c; // move to next element of local matrix
}
++t_row; // move to next base function on row
}
++t_w; // move to next integrate weight
}
PetscFunctionReturn(0);
}
/**
* \brief integrate local constrains matrix
*/
inline PetscErrorCode aSsemble(
DataForcesAndSourcesCore::EntData &row_data,DataForcesAndSourcesCore::EntData &col_data
) {
PetscFunctionBegin;
// get indices on row
const int* row_indices = &*row_data.getIndices().data().begin();
// get indices on column
const int* col_indices = &*col_data.getIndices().data().begin();
Mat B = getFEMethod()->ksp_B!=PETSC_NULL? getFEMethod()->ksp_B : getFEMethod()->snes_B;
// assemble local matrix
ierr = MatSetValues(
B,nbRows,row_indices,nbCols,col_indices,&*locMat.data().begin(),ADD_VALUES
// cerr << locMat << endl;
if(assembleTraspose) {
// assmble transpose of local matrix
locMat = trans(locMat);
ierr = MatSetValues(
B, nbCols,col_indices,nbRows,row_indices,&*locMat.data().begin(),ADD_VALUES
}
PetscFunctionReturn(0);
}
};
/**
* \brief Assemble constrains vector
*
* \f[
* \mathbf{g} = \int_{\partial\Omega} \boldsymbol\psi \overline{u} \textrm{d}\partial\Omega
* \f]
*
*/
struct Op_g: public OpBaseRhs<FaceElementForcesAndSourcesCore::UserDataOperator> {
typedef boost::function<double (const double,const double,const double)> FVal;
Op_g(FVal f_value,const string field_name = "L",const double beta = 1):
OpBaseRhs<FaceElementForcesAndSourcesCore::UserDataOperator>(field_name),
fValue(f_value),
bEta(beta) {
}
protected:
FTensor::Number<0> NX; ///< x-direction index
FTensor::Number<1> NY; ///< y-direction index
FTensor::Number<2> NZ; ///< z-direction index
FVal fValue; ///< Function pointer evaluating values of "U" at the boundary
VectorDouble locVec;
const double bEta;
/**
* \brief Integrate local constrains vector
*/
PetscErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &data) {
PetscFunctionBegin;
// set size to local vector
locVec.resize(nbRows,false);
// clear loacl vector
locVec.clear();
// get face area
const double area = getArea()*bEta;
// get integration wiegth
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
// get base function
FTensor::Tensor0<double*> t_l = data.getFTensor0N();
// get coordinates at integration point
FTensor::Tensor1<double*,3> t_coords = getTensor1CoordsAtGaussPts();
// make loop over integration points
for(int gg = 0;gg!=nbIntegrationPts;gg++) {
// evalue function on boundary and scale it by area and integration weight
double alpha = area*t_w*fValue(t_coords(NX),t_coords(NY),t_coords(NZ));
// get element of vector
FTensor::Tensor0<double*> t_a(&*locVec.data().begin());
//
for(int rr = 0;rr!=nbRows;rr++) {
t_a += alpha*t_l;
++t_a;
++t_l;
}
++t_w;
++t_coords;
}
PetscFunctionReturn(0);
}
/**
* \brief assemble constrains vectors
*/
PetscErrorCode aSsemble(DataForcesAndSourcesCore::EntData &data) {
PetscFunctionBegin;
const int* indices = &*data.getIndices().data().begin();
const double* vals = &*locVec.data().begin();
Vec f = getFEMethod()->ksp_f!=PETSC_NULL ? getFEMethod()->ksp_f : getFEMethod()->snes_f;
ierr = VecSetValues(f,nbRows,indices,&*vals,ADD_VALUES); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
};
/**
* \brief Evaluate error
*/
struct OpError: public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
typedef boost::function<double (const double,const double,const double)> UVal;
typedef boost::function<FTensor::Tensor1<double,3> (const double,const double,const double)> GVal;
OpError(
UVal u_value,
GVal g_value,
boost::shared_ptr<VectorDouble>& field_vals,
boost::shared_ptr<MatrixDouble>& grad_vals,
Vec global_error
):
OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator>("ERROR"),
globalError(global_error),
uValue(u_value),
gValue(g_value),
fieldVals(field_vals),
gradVals(grad_vals) {
}
PetscErrorCode doWork(
int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
) {
PetscFunctionBegin;
nbRows = row_data.getFieldData().size();
if(!nbRows) PetscFunctionReturn(0);
nbIntegrationPts = getGaussPts().size2();
ierr = iNtegrate(row_data); CHKERRQ(ierr);
ierr = aSsemble(row_data); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
private:
Vec globalError; ///< ghost vector with global (integrated over volume) error
UVal uValue; ///< function with exact solution
GVal gValue; ///< function with exact solution for gradient
boost::shared_ptr<VectorDouble> fieldVals;
boost::shared_ptr<MatrixDouble> gradVals;
/**
* \brief Integrate error
*/
PetscErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &data) {
PetscFunctionBegin;
// clear field dofs
data.getFieldData().clear();
// get volume of element
const double vol = getVolume();
// get integration weight
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
// get solution at integration point
// get solution at integration point
FTensor::Tensor1<double*,3> t_grad = getTensor1FormData<3>(*gradVals);
// get coordinates at integration point
FTensor::Tensor1<double*,3> t_coords = getTensor1CoordsAtGaussPts();
// keep exact gradient and error or gradient
FTensor::Tensor1<double,3> t_exact_grad,t_error_grad;
// integrate over
for(int gg = 0;gg!=nbIntegrationPts;gg++) {
double alpha = vol*t_w;
// evalue exact value
double exact_u = uValue(t_coords(NX),t_coords(NY),t_coords(NZ));
// evalue exact hradient
t_exact_grad = gValue(t_coords(NX),t_coords(NY),t_coords(NZ));
// calculate gradient errro
t_error_grad(i) = t_grad(i)-t_exact_grad(i);
// error
double error = pow(t_u-exact_u,2)+t_error_grad(i)*t_error_grad(i);
// iterate over base functions
data.getFieldData()[0] += alpha*error;
++t_w; // move to next integration point
++t_u; // next value of function at integration point
++t_grad; // next gradient at integration point
++t_coords; // next coordinate at integration point
}
PetscFunctionReturn(0);
}
/**
* \brief Assemble error
*/
PetscErrorCode aSsemble(DataForcesAndSourcesCore::EntData &data) {
PetscFunctionBegin;
// set error on mesh
data.getFieldDofs()[0]->getFieldData() = sqrt(data.getFieldData()[0]);
// assemble vector to global error
ierr = VecSetValue(globalError,0,data.getFieldData()[0],ADD_VALUES); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
};
struct OpKt: public OpK {
OpKt(
boost::function<double (const double)> a,
boost::function<double (const double)> diff_a,
boost::shared_ptr<VectorDouble>& field_vals,
boost::shared_ptr<MatrixDouble>& grad_vals
):
OpK(false),
A(a),
diffA(diff_a),
fieldVals(field_vals),
gradVals(grad_vals) {
}
protected:
/**
* \brief Integrate grad-grad operator
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
inline PetscErrorCode iNtegrate(
DataForcesAndSourcesCore::EntData &row_data,DataForcesAndSourcesCore::EntData &col_data
) {
PetscFunctionBegin;
// set size of local entity bock
locMat.resize(nbRows,nbCols,false);
// clear matrix
locMat.clear();
// get element volume
double vol = getVolume();
// get integration weights
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
// get solution at integration point
// get solution at integration point
FTensor::Tensor1<double*,3> t_grad = getTensor1FormData<3>(*gradVals);
// get base function gradient on rows
FTensor::Tensor1<double*,3> t_row_grad = row_data.getFTensor1DiffN<3>();
// loop over integration points
for(int gg = 0;gg!=nbIntegrationPts;gg++) {
// take into account Jacobian
const double alpha = t_w*vol;
const double beta = alpha*A(t_u);
t_gamma(i) = (alpha*diffA(t_u))*t_grad(i);
// take fist element to local matrix
FTensor::Tensor0<double*> a(&*locMat.data().begin());
// loop over rows base functions
for(int rr = 0;rr!=nbRows;rr++) {
// get column base function
FTensor::Tensor0<double*> t_col = col_data.getFTensor0N(gg,0);
// get column base functions gradient at gauss point gg
FTensor::Tensor1<double*,3> t_col_grad = col_data.getFTensor1DiffN<3>(gg,0);
// loop over columns
for(int cc = 0;cc!=nbCols;cc++) {
// calculate element of local matrix
a += (t_row_grad(i)*beta)*t_col_grad(i)+t_row_grad(i)*(t_gamma(i)*t_col);
++t_col; // move to next base function
++t_col_grad; // move to another gradient of base function on column
++a; // move to another element of local matrix in column
}
++t_row_grad; // move to another element of gradient of base function on row
}
++t_w; // move to another integration weight
++t_u; // move to next value at integration point
++t_grad; // move to next gradient value
}
PetscFunctionReturn(0);
}
boost::function<double (const double)> A;
boost::function<double (const double)> diffA;
boost::shared_ptr<VectorDouble> fieldVals;
boost::shared_ptr<MatrixDouble> gradVals;
};
struct OpResF_Domain: public OpF {
OpResF_Domain(
FSource f_source,
boost::function<double (const double)> a,
boost::shared_ptr<VectorDouble>& field_vals,
boost::shared_ptr<MatrixDouble>& grad_vals
):
OpF(f_source),
A(a),
fieldVals(field_vals),
gradVals(grad_vals) {
}
protected:
/**
* \brief Integrate local entity vector
* @param data entity data on element row
* @return error code
*/
PetscErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &data) {
PetscFunctionBegin;
// set size of local vector
locVec.resize(nbRows,false);
// clear local entity vector
locVec.clear();
// get finite element volume
double vol = getVolume();
// get integration weights
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
// get solution at integration point
// get solution at integration point
FTensor::Tensor1<double*,3> t_grad = getTensor1FormData<3>(*gradVals);
// get base functions on entity
FTensor::Tensor0<double*> t_v = data.getFTensor0N();
// get base function gradient on rows
FTensor::Tensor1<double*,3> t_v_grad = data.getFTensor1DiffN<3>();
// get coordinates at integration points
FTensor::Tensor1<double*,3> t_coords = getTensor1CoordsAtGaussPts();
// loop over all integration points
for(int gg = 0;gg!=nbIntegrationPts;gg++) {
// evaluate constant term
const double alpha = vol*t_w;
const double source_term = alpha*fSource(t_coords(NX),t_coords(NY),t_coords(NZ));
grad_term(i) = (alpha*A(t_u))*t_grad(i);
// get element of local vector
FTensor::Tensor0<double*> t_a(&*locVec.data().begin());
// loop over base functions
for(int rr = 0;rr!=nbRows;rr++) {
// add to local vector source term
t_a += t_v_grad(i)*grad_term(i)+t_v*source_term;
++t_a; // move to next element of local vector
++t_v; // move to next base function
++t_v_grad; // move to next gradient of base function
}
++t_w; // move to next integration weights
++t_u; // move to next value
++t_grad; // move to next gradient value
++t_coords; // move to next physical coordinates at integration point
}
PetscFunctionReturn(0);
}
boost::function<double (const double)> A;
boost::shared_ptr<VectorDouble> fieldVals;
boost::shared_ptr<MatrixDouble> gradVals;
};
struct OpRes_g: public Op_g {
OpRes_g(
FVal f_value,
boost::shared_ptr<VectorDouble>& field_vals
):
Op_g(f_value,"L",1),
fieldVals(field_vals) {
}
protected:
boost::shared_ptr<VectorDouble> fieldVals;
/**
* \brief Integrate local constrains vector
*/
PetscErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &data) {
PetscFunctionBegin;
// set size to local vector
locVec.resize(nbRows,false);
// clear loacl vector
locVec.clear();
// get face area
const double area = getArea()*bEta;
// get integration wiegth
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
// get base function
FTensor::Tensor0<double*> t_l = data.getFTensor0N();
// get solution at integration point
// get coordinates at integration point
FTensor::Tensor1<double*,3> t_coords = getTensor1CoordsAtGaussPts();
// make loop over integration points
for(int gg = 0;gg!=nbIntegrationPts;gg++) {
// evalue function on boundary and scale it by area and integration weight
double alpha = area*t_w;
// get element of vector
FTensor::Tensor0<double*> t_a(&*locVec.data().begin());
for(int rr = 0;rr!=nbRows;rr++) {
t_a += alpha*t_l*(t_u-fValue(t_coords(NX),t_coords(NY),t_coords(NZ)));
++t_a;
++t_l;
}
++t_w;
++t_u;
++t_coords;
}
PetscFunctionReturn(0);
}
};
struct OpResF_Boundary: public Op_g {
OpResF_Boundary(
boost::shared_ptr<VectorDouble>& lambda_vals
):
Op_g(FVal(),"U",1),
lambdaVals(lambda_vals) {
}
protected:
boost::shared_ptr<VectorDouble> lambdaVals;
/**
* \brief Integrate local constrains vector
*/
PetscErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &data) {
PetscFunctionBegin;
// set size to local vector
locVec.resize(nbRows,false);
// clear loacl vector
locVec.clear();
// get face area
const double area = getArea()*bEta;
// get integration wiegth
FTensor::Tensor0<double*> t_w = getFTensor0IntegrationWeight();
// get base function
FTensor::Tensor0<double*> t_u = data.getFTensor0N();
// get solution at integration point
// make loop over integration points
for(int gg = 0;gg!=nbIntegrationPts;gg++) {
// evalue function on boundary and scale it by area and integration weight
double alpha = area*t_w;
// get element of vector
FTensor::Tensor0<double*> t_a(&*locVec.data().begin());
for(int rr = 0;rr!=nbRows;rr++) {
t_a += alpha*t_u*t_lambda;
++t_a;
++t_u;
}
++t_w;
++t_lambda;
}
PetscFunctionReturn(0);
}
};
/**
* \brief Set integration rule to volume elements
*
* This rule is used to integrate \f$\nabla v \cdot \nabla u\f$, thus
* if approximation field and testing field are polynomial order "p",
* then rule for exact integration is 2*(p-1).
*
* Integration rule is order of polynomial which is calculated exactly. Finite element
* selects integration method based on return of this function.
*
*/
struct VolRule {
int operator()(int,int,int p) const {
return 2*(p-1);
}
};
/**
* \brief Set integration rule to boundary elements
*
* This is uses to integrate values on the face. Is used to integrate
* \f$(\mathbf{n} \cdot \lambda) u\f$, where Lagrange multiplayer
* is order "p_row" and approximate function is order "p_col".
*
* Integration rule is order of polynomial which is calculated exactly. Finite element
* selects integration method based on return of this function.
*
*/
struct FaceRule {
int operator()(int p_row,int p_col,int p_data) const {
return 2*p_data+1;
}
};
/**
* \brief Create finite elements instances
*
* Create finite element instances and add operators to finite elements.
*
*/
struct CreateFiniteElements {
CreateFiniteElements(MoFEM::Interface &m_field):
mField(m_field) {
}
/**
* \brief Create finite element to calculate matrix and vectors
*/
PetscErrorCode createFEToAssmbleMatrceAndVector(
boost::function<double (const double,const double,const double)> f_u,
boost::function<double (const double,const double,const double)> f_source,
boost::shared_ptr<ForcesAndSourcesCore>& domain_lhs_fe,
boost::shared_ptr<ForcesAndSourcesCore>& boundary_lhs_fe,
boost::shared_ptr<ForcesAndSourcesCore>& domain_rhs_fe,
boost::shared_ptr<ForcesAndSourcesCore>& boundary_rhs_fe,
bool trans = true
) const {
PetscFunctionBegin;
// Create elements element instances
domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(new VolumeElementForcesAndSourcesCore(mField));
boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(new FaceElementForcesAndSourcesCore(mField));
domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(new VolumeElementForcesAndSourcesCore(mField));
boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(new FaceElementForcesAndSourcesCore(mField));
// Set integration rule to elements instances
domain_lhs_fe->getRuleHook = VolRule();
domain_rhs_fe->getRuleHook = VolRule();
boundary_lhs_fe->getRuleHook = FaceRule();
boundary_rhs_fe->getRuleHook = FaceRule();
// Add operators to element instances
// Add operator grad-grad for calculaye matrix
domain_lhs_fe->getOpPtrVector().push_back(new OpK());
// Add operator to calculate source terms
domain_rhs_fe->getOpPtrVector().push_back(new OpF(f_source));
// Add operator calculating constrains matrix
boundary_lhs_fe->getOpPtrVector().push_back(new OpC(trans));
// Add operator calculating constrains vector
boundary_rhs_fe->getOpPtrVector().push_back(new Op_g(f_u));
PetscFunctionReturn(0);
}
/**
* \brief Create finite element to calculate error
*/
PetscErrorCode createFEToEvaluateError(
boost::function<double (const double,const double,const double)> f_u,
boost::function<FTensor::Tensor1<double,3> (const double,const double,const double)> g_u,
Vec global_error,
boost::shared_ptr<ForcesAndSourcesCore>& domain_error
) const {
PetscFunctionBegin;
// Create finite element instance to calculaye error
domain_error = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField)
);
domain_error->getRuleHook = VolRule();
// Set integration rule
// Crate shared vector storing values of field "u" on integration points on element. element
// is local and is used to exchange data between operators.
boost::shared_ptr<VectorDouble> values_at_integation_ptr = boost::make_shared<VectorDouble>();
// Storing gradients of field
boost::shared_ptr<MatrixDouble> grad_at_integation_ptr = boost::make_shared<MatrixDouble>();
// Add default operator to calculate field values at integration points
domain_error->getOpPtrVector().push_back(new OpCalculateScalarFieldValues("U",values_at_integation_ptr));
// Add default operator to calculate field gradient at integration points
domain_error->getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>("U",grad_at_integation_ptr));
// Add operator to integrate error element by element.
domain_error->getOpPtrVector().push_back(
new OpError(f_u,g_u,values_at_integation_ptr,grad_at_integation_ptr,global_error)
);
PetscFunctionReturn(0);
}
/**
* \brief Create finite element to post-process results
*/
PetscErrorCode creatFEToPostProcessResults(
boost::shared_ptr<ForcesAndSourcesCore>& post_proc_volume
) const {
PetscFunctionBegin;
// Note that user can stack together arbitrary number of operators to compose
// complex PDEs.
// Post-process results. This is standard element, with functionality
// enabling refining mesh for post-processing. In addition in PostProcOnRefMesh.hpp
// are implanted set of users operators to post-processing fields. Here
// using simplified mechanism for post-processing finite element, we
// add operators to save data from field on mesh tags for ParaView
// visualization.
post_proc_volume = boost::shared_ptr<ForcesAndSourcesCore>(new PostProcVolumeOnRefinedMesh(mField));
// Add operators to the elements, starting with some generic
ierr = boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(post_proc_volume)->
generateReferenceElementMesh(); CHKERRQ(ierr);
ierr = boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(post_proc_volume)->
addFieldValuesPostProc("U"); CHKERRQ(ierr);
ierr = boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(post_proc_volume)->
addFieldValuesPostProc("ERROR"); CHKERRQ(ierr);
ierr = boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(post_proc_volume)->
addFieldValuesGradientPostProc("U"); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
/**
* \brief Create finite element to calculate matrix and vectors
*/
PetscErrorCode createFEToAssmbleMatrceAndVectorForNonlinearProblem(
boost::function<double (const double,const double,const double)> f_u,
boost::function<double (const double,const double,const double)> f_source,
boost::function<double (const double)> a,
boost::function<double (const double)> diff_a,
boost::shared_ptr<ForcesAndSourcesCore>& domain_lhs_fe,
boost::shared_ptr<ForcesAndSourcesCore>& boundary_lhs_fe,
boost::shared_ptr<ForcesAndSourcesCore>& domain_rhs_fe,
boost::shared_ptr<ForcesAndSourcesCore>& boundary_rhs_fe,
ForcesAndSourcesCore::RuleHookFun vol_rule,
ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule(),
bool trans = true
) const {
PetscFunctionBegin;
// Create elements element instances
domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(new VolumeElementForcesAndSourcesCore(mField));
boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(new FaceElementForcesAndSourcesCore(mField));
domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(new VolumeElementForcesAndSourcesCore(mField));
boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(new FaceElementForcesAndSourcesCore(mField));
// Set integration rule to elements instances
domain_lhs_fe->getRuleHook = vol_rule;
domain_rhs_fe->getRuleHook = vol_rule;
boundary_lhs_fe->getRuleHook = face_rule;
boundary_rhs_fe->getRuleHook = face_rule;
// Set integration rule
// Crate shared vector storing values of field "u" on integration points on element. element
// is local and is used to exchange data between operators.
boost::shared_ptr<VectorDouble> values_at_integation_ptr = boost::make_shared<VectorDouble>();
// Storing gradients of field
boost::shared_ptr<MatrixDouble> grad_at_integation_ptr = boost::make_shared<MatrixDouble>();
// multipliers values
boost::shared_ptr<VectorDouble> multiplier_at_integation_ptr = boost::make_shared<VectorDouble>();
// Add operators to element instances
// Add default operator to calculate field values at integration points
domain_lhs_fe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues("U",values_at_integation_ptr));
// Add default operator to calculate field gradient at integration points
domain_lhs_fe->getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>("U",grad_at_integation_ptr));
// Add operator grad-(1+u^2)grad for calculate matrix
domain_lhs_fe->getOpPtrVector().push_back(new OpKt(a,diff_a,values_at_integation_ptr,grad_at_integation_ptr));
// Add default operator to calculate field values at integration points
domain_rhs_fe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues("U",values_at_integation_ptr));
// Add default operator to calculate field gradient at integration points
domain_rhs_fe->getOpPtrVector().push_back(new OpCalculateScalarFieldGradient<3>("U",grad_at_integation_ptr));
// Add operator to calculate source terms
domain_rhs_fe->getOpPtrVector().push_back(new OpResF_Domain(f_source,a,values_at_integation_ptr,grad_at_integation_ptr));
// Add operator calculating constrains matrix
boundary_lhs_fe->getOpPtrVector().push_back(new OpC(trans));
// Add default operator to calculate field values at integration points
boundary_rhs_fe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues("U",values_at_integation_ptr));
// Add defualy operator to calculate values of Lagrange multipliers
boundary_rhs_fe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues("L",multiplier_at_integation_ptr));
// Add operator calculating constrains vector
boundary_rhs_fe->getOpPtrVector().push_back(new OpRes_g(f_u,values_at_integation_ptr));
boundary_rhs_fe->getOpPtrVector().push_back(new OpResF_Boundary(multiplier_at_integation_ptr));
PetscFunctionReturn(0);
}
private:
};
}
#endif //__POISSONOPERATORS_HPP__