Implementation of operators, problem and finite elements for unsaturated flow. More...

#include "tutorials/cor-0to1/src/UnsaturatedFlow.hpp"

Inheritance diagram for MixTransport::UnsaturatedFlowElement:

Collaboration diagram for MixTransport::UnsaturatedFlowElement:

Classes
struct	BcData
	Class storing information about boundary condition. More...

struct	FaceRule
	Set integration rule to boundary elements. More...

struct	MonitorPostProc

struct	OpDivTauU_HdivL2

struct	OpEvaluateInitiallHead

struct	OpIntegrateFluxes

struct	OpPostProcMaterial

struct	OpResidualFlux
	Assemble flux residual. More...

struct	OpResidualMass

struct	OpRhsBcOnValues
	Evaluate boundary condition at the boundary. More...

struct	OpTauDotSigma_HdivHdiv

struct	OpVDivSigma_L2Hdiv

struct	OpVU_L2L2

struct	postProcessVol
	Post proces method for volume element Assemble vectors and matrices and apply essential boundary conditions. More...

struct	preProcessVol
	Pre-peprocessing Set head pressure rate and get inital essential boundary conditions. More...

struct	VolRule
	Set integration rule to volume elements. More...

Public Types
typedef std::map< int, boost::shared_ptr< GenericMaterial > >	MaterialsDoubleMap

typedef map< int, boost::shared_ptr< BcData > >	BcMap

Public Member Functions
	UnsaturatedFlowElement (MoFEM::Interface &m_field)

	~UnsaturatedFlowElement ()

virtual MoFEMErrorCode	getMaterial (const EntityHandle ent, int &block_id) const
	For given element handle get material block Id.

MoFEMErrorCode	getBcOnValues (const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
	Get value on boundary.

MoFEMErrorCode	getBcOnFluxes (const EntityHandle ent, const double x, const double y, const double z, double &flux)
	essential (Neumann) boundary condition (set fluxes)

MoFEMErrorCode	addFields (const std::string &values, const std::string &fluxes, const int order)
	add fields

MoFEMErrorCode	addFiniteElements (const std::string &fluxes_name, const std::string &values_name)
	add finite elements

MoFEMErrorCode	buildProblem (Range zero_flux_ents, BitRefLevel ref_level=BitRefLevel().set(0))
	Build problem.

MoFEMErrorCode	setFiniteElements (ForcesAndSourcesCore::RuleHookFun vol_rule=VolRule(), ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule())
	Create finite element instances.

MoFEMErrorCode	createMatrices ()
	Create vectors and matrices.

MoFEMErrorCode	destroyMatrices ()
	Delete matrices and vector when no longer needed.

MoFEMErrorCode	calculateEssentialBc ()
	Calculate boundary conditions for fluxes.

MoFEMErrorCode	calculateInitialPc ()
	Calculate inital pressure head distribution.

MoFEMErrorCode	solveProblem (bool set_initial_pc=true)
	solve problem

Public Member Functions inherited from MixTransport::MixTransportElement
	MixTransportElement (MoFEM::Interface &m_field)
	construction of this data structure

virtual	~MixTransportElement ()
	destructor

MoFEMErrorCode	getDirichletBCIndices (IS *is)
	get dof indices where essential boundary conditions are applied

virtual MoFEMErrorCode	getSource (const EntityHandle ent, const double x, const double y, const double z, double &flux)
	set source term

virtual MoFEMErrorCode	getResistivity (const EntityHandle ent, const double x, const double y, const double z, MatrixDouble3by3 &inv_k)
	natural (Dirichlet) boundary conditions (set values)

MoFEMErrorCode	addFields (const std::string &values, const std::string &fluxes, const int order)
	Add fields to database.

MoFEMErrorCode	addFiniteElements (const std::string &fluxes_name, const std::string &values_name)
	add finite elements

MoFEMErrorCode	buildProblem (BitRefLevel &ref_level)
	Build problem.

MoFEMErrorCode	postProc (const string out_file)
	Post process results.

MoFEMErrorCode	createMatrices ()
	create matrices

MoFEMErrorCode	solveLinearProblem ()
	solve problem

MoFEMErrorCode	calculateResidual ()
	calculate residual

MoFEMErrorCode	evaluateError ()
	Calculate error on elements.

MoFEMErrorCode	destroyMatrices ()
	destroy matrices

Public Attributes
DM	dM
	Discrete manager for unsaturated flow problem.

MaterialsDoubleMap	dMatMap
	materials database

BcMap	bcValueMap
	Store boundary condition for head capillary pressure.

EntityHandle	lastEvalBcValEnt

int	lastEvalBcBlockValId

BcMap	bcFluxMap

EntityHandle	lastEvalBcFluxEnt

int	lastEvalBcBlockFluxId

boost::shared_ptr< ForcesAndSourcesCore >	feFaceBc
	Elemnet to calculate essential bc.

boost::shared_ptr< ForcesAndSourcesCore >	feFaceRhs
	Face element apply natural bc.

boost::shared_ptr< ForcesAndSourcesCore >	feVolInitialPc
	Calculate inital boundary conditions.

boost::shared_ptr< ForcesAndSourcesCore >	feVolRhs
	Assemble residual vector.

boost::shared_ptr< ForcesAndSourcesCore >	feVolLhs
	Assemble tangent matrix.

boost::shared_ptr< MethodForForceScaling >	scaleMethodFlux
	Method scaling fluxes.

boost::shared_ptr< MethodForForceScaling >	scaleMethodValue
	Method scaling values.

boost::shared_ptr< FEMethod >	tsMonitor

boost::shared_ptr< VectorDouble >	headRateAtGaussPts
	Vector keeps head rate.

std::vector< int >	bcVecIds

VectorDouble	bcVecVals

VectorDouble	vecValsOnBc

Vec	D1
	Vector with inital head capillary pressure.

Vec	ghostFlux
	Ghost Vector of integrated fluxes.

Public Attributes inherited from MixTransport::MixTransportElement
MoFEM::Interface &	mField

MyVolumeFE	feVol
	Instance of volume element.

MyTriFE	feTri
	Instance of surface element.

VectorDouble	valuesAtGaussPts
	values at integration points on element

MatrixDouble	valuesGradientAtGaussPts
	gradients at integration points on element

VectorDouble	divergenceAtGaussPts
	divergence at integration points on element

MatrixDouble	fluxesAtGaussPts
	fluxes at integration points on element

set< int >	bcIndices

std::map< int, BlockData >	setOfBlocks
	maps block set id with appropriate BlockData

Vec	D

Vec	D0

Vec	F

Mat	Aij

map< double, EntityHandle >	errorMap

double	sumErrorFlux

double	sumErrorDiv

double	sumErrorJump

Detailed Description

Implementation of operators, problem and finite elements for unsaturated flow.

Examples: unsaturated_transport.cpp.

Definition at line 102 of file UnsaturatedFlow.hpp.

Member Typedef Documentation

◆ BcMap

map<int, boost::shared_ptr<BcData> > MixTransport::UnsaturatedFlowElement::BcMap

Examples: UnsaturatedFlow.hpp.

Definition at line 152 of file UnsaturatedFlow.hpp.

◆ MaterialsDoubleMap

std::map<int, boost::shared_ptr<GenericMaterial> > MixTransport::UnsaturatedFlowElement::MaterialsDoubleMap

Examples: UnsaturatedFlow.hpp.

Definition at line 118 of file UnsaturatedFlow.hpp.

Constructor & Destructor Documentation

◆ UnsaturatedFlowElement()

MixTransport::UnsaturatedFlowElement::UnsaturatedFlowElement ( MoFEM::Interface & m_field )

inline

Examples: UnsaturatedFlow.hpp.

Definition at line 106 of file UnsaturatedFlow.hpp.

      : MixTransportElement(m_field), dM(PETSC_NULLPTR), lastEvalBcValEnt(0),
        lastEvalBcBlockValId(-1), lastEvalBcFluxEnt(0),
        lastEvalBcBlockFluxId(-1) {}

◆ ~UnsaturatedFlowElement()

MixTransport::UnsaturatedFlowElement::~UnsaturatedFlowElement ( )

inline

Examples: UnsaturatedFlow.hpp.

Definition at line 111 of file UnsaturatedFlow.hpp.

                            {
    if (dM != PETSC_NULLPTR) {
      CHKERR DMDestroy(&dM);
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
    }
  }

Member Function Documentation

◆ addFields()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::addFields	(	const std::string &	values,
		const std::string &	fluxes,
		const int	order )

inline

add fields

Examples: UnsaturatedFlow.hpp, and unsaturated_transport.cpp.

Definition at line 1136 of file UnsaturatedFlow.hpp.

                                            {
    MoFEMFunctionBegin;
    // Fields
    CHKERR mField.add_field(fluxes, HDIV, DEMKOWICZ_JACOBI_BASE, 1);
    CHKERR mField.add_field(values, L2, AINSWORTH_LEGENDRE_BASE, 1);
    CHKERR mField.add_field(values + "_t", L2, AINSWORTH_LEGENDRE_BASE, 1);
    // CHKERR mField.add_field(fluxes+"_residual",L2,AINSWORTH_LEGENDRE_BASE,1);
 
    // meshset consisting all entities in mesh
    EntityHandle root_set = mField.get_moab().get_root_set();
    // add entities to field
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, 4, "SOIL") != 0)
        continue;
      CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
                                              MBTET, fluxes);
      CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
                                              MBTET, values);
      CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
                                              MBTET, values + "_t");
      // CHKERR mField.add_ents_to_field_by_type(
      //   dMatMap[it->getMeshsetId()]->tEts,MBTET,fluxes+"_residual"
      // );
    }
 
    CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
    CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
    CHKERR mField.set_field_order(root_set, MBTET, values, order);
    CHKERR mField.set_field_order(root_set, MBTET, values + "_t", order);
    // CHKERR mField.set_field_order(root_set,MBTET,fluxes+"_residual",order);
    MoFEMFunctionReturn(0);
  }

◆ addFiniteElements()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::addFiniteElements	(	const std::string &	fluxes_name,
		const std::string &	values_name )

inline

add finite elements

Examples: UnsaturatedFlow.hpp, and unsaturated_transport.cpp.

Definition at line 1172 of file UnsaturatedFlow.hpp.

                                                                 {
    MoFEMFunctionBegin;
 
    // Define element "MIX". Note that this element will work with fluxes_name
    // and values_name. This reflect bilinear form for the problem
    CHKERR mField.add_finite_element("MIX", MF_ZERO);
    CHKERR mField.modify_finite_element_add_field_row("MIX", fluxes_name);
    CHKERR mField.modify_finite_element_add_field_col("MIX", fluxes_name);
    CHKERR mField.modify_finite_element_add_field_row("MIX", values_name);
    CHKERR mField.modify_finite_element_add_field_col("MIX", values_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX", fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX", values_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX",
                                                       values_name + "_t");
    // CHKERR
    // mField.modify_finite_element_add_field_data("MIX",fluxes_name+"_residual");
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, 4, "SOIL") != 0)
        continue;
      CHKERR mField.add_ents_to_finite_element_by_type(
          dMatMap[it->getMeshsetId()]->tEts, MBTET, "MIX");
    }
 
    // Define element to integrate natural boundary conditions, i.e. set values.
    CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
    CHKERR mField.modify_finite_element_add_field_row("MIX_BCVALUE",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_col("MIX_BCVALUE",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX_BCVALUE",
                                                       fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX_BCVALUE",
                                                       values_name);
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, 4, "HEAD") != 0)
        continue;
      CHKERR mField.add_ents_to_finite_element_by_type(
          bcValueMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCVALUE");
    }
 
    // Define element to apply essential boundary conditions.
    CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
    CHKERR mField.modify_finite_element_add_field_row("MIX_BCFLUX",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_col("MIX_BCFLUX",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX_BCFLUX",
                                                       fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX_BCFLUX",
                                                       values_name);
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, 4, "FLUX") != 0)
        continue;
      CHKERR mField.add_ents_to_finite_element_by_type(
          bcFluxMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCFLUX");
    }
 
    MoFEMFunctionReturn(0);
  }

◆ buildProblem()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::buildProblem	(	Range	zero_flux_ents,
		BitRefLevel	ref_level = BitRefLevel().set(0) )

inline

Build problem.

Parameters

ref_level mesh refinement on which mesh problem you like to built.

Returns: error code

Examples: UnsaturatedFlow.hpp, and unsaturated_transport.cpp.

Definition at line 1241 of file UnsaturatedFlow.hpp.

                                                                            {
    MoFEMFunctionBegin;
 
    // Build fields
    CHKERR mField.build_fields();
    // Build finite elements
    CHKERR mField.build_finite_elements("MIX");
    CHKERR mField.build_finite_elements("MIX_BCFLUX");
    CHKERR mField.build_finite_elements("MIX_BCVALUE");
    // Build adjacencies of degrees of freedom and elements
    CHKERR mField.build_adjacencies(ref_level);
 
    //  create DM instance
    CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
    // setting that DM is type of DMMOFEM, i.e. MOFEM implementation manages DM
    CHKERR DMSetType(dM, "DMMOFEM");
    // mesh is portioned, each process keeps only part of problem
    CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
    // creates problem in DM
    CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
    // discretised problem creates square matrix (that makes some optimizations)
    CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
    // set DM options from command line
    CHKERR DMSetFromOptions(dM);
    // add finite elements
    CHKERR DMMoFEMAddElement(dM, "MIX");
    CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
    CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
    // constructor data structures
    CHKERR DMSetUp(dM);
 
    // remove zero flux dofs
    CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
        "MIX", "FLUXES", zero_flux_ents);
 
    PetscSection section;
    CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
    CHKERR DMSetSection(dM, section);
    // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
    CHKERR PetscSectionDestroy(&section);
 
    MoFEMFunctionReturn(0);
  }

◆ calculateEssentialBc()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::calculateEssentialBc ( )

inline

Calculate boundary conditions for fluxes.

Returns: Error code

Examples: unsaturated_transport.cpp.

Definition at line 1649 of file UnsaturatedFlow.hpp.

                                        {
    MoFEMFunctionBegin;
    // clear vectors
    CHKERR VecZeroEntries(D0);
    CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
    // clear essential bc indices, it could have dofs from other mesh refinement
    bcIndices.clear();
    // set operator to calculate essential boundary conditions
    CHKERR DMoFEMLoopFiniteElements(dM, "MIX_BCFLUX", feFaceBc);
    CHKERR VecAssemblyBegin(D0);
    CHKERR VecAssemblyEnd(D0);
    CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
    double norm2D0;
    CHKERR VecNorm(D0, NORM_2, &norm2D0);
    // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
    PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
    MoFEMFunctionReturn(0);
  }

◆ calculateInitialPc()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::calculateInitialPc ( )

inline

Calculate inital pressure head distribution.

Returns: Error code

Examples: unsaturated_transport.cpp.

Definition at line 1674 of file UnsaturatedFlow.hpp.

                                      {
    MoFEMFunctionBegin;
    // clear vectors
    CHKERR VecZeroEntries(D1);
    CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
    // Calculate initial pressure head on each element
    CHKERR DMoFEMLoopFiniteElements(dM, "MIX", feVolInitialPc);
    // Assemble vector
    CHKERR VecAssemblyBegin(D1);
    CHKERR VecAssemblyEnd(D1);
    CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
    // Calculate norm
    double norm2D1;
    CHKERR VecNorm(D1, NORM_2, &norm2D1);
    // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
    PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
    MoFEMFunctionReturn(0);
  }

◆ createMatrices()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::createMatrices ( )

inline

Create vectors and matrices.

Returns: Error code

Examples: unsaturated_transport.cpp.

Definition at line 1615 of file UnsaturatedFlow.hpp.

                                  {
    MoFEMFunctionBegin;
    CHKERR DMCreateMatrix(dM, &Aij);
    CHKERR DMCreateGlobalVector(dM, &D0);
    CHKERR VecDuplicate(D0, &D1);
    CHKERR VecDuplicate(D0, &D);
    CHKERR VecDuplicate(D0, &F);
    int ghosts[] = {0};
    int nb_locals = mField.get_comm_rank() == 0 ? 1 : 0;
    int nb_ghosts = mField.get_comm_rank() > 0 ? 1 : 0;
    CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
                          &ghostFlux);
    MoFEMFunctionReturn(0);
  }

◆ destroyMatrices()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::destroyMatrices ( )

inline

Delete matrices and vector when no longer needed.

Returns: error code

Examples: unsaturated_transport.cpp.

Definition at line 1634 of file UnsaturatedFlow.hpp.

                                   {
    MoFEMFunctionBegin;
    CHKERR MatDestroy(&Aij);
    CHKERR VecDestroy(&D);
    CHKERR VecDestroy(&D0);
    CHKERR VecDestroy(&D1);
    CHKERR VecDestroy(&F);
    CHKERR VecDestroy(&ghostFlux);
    MoFEMFunctionReturn(0);
  }

◆ getBcOnFluxes()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::getBcOnFluxes	(	const EntityHandle	ent,
		const double	x,
		const double	y,
		const double	z,
		double &	flux )

inlinevirtual

essential (Neumann) boundary condition (set fluxes)

Parameters

ent	handle to finite element entity
x	coord
y	coord
z	coord
flux	reference to flux which is set by function

Returns: [description]

Reimplemented from MixTransport::MixTransportElement.

Examples: UnsaturatedFlow.hpp.

Definition at line 210 of file UnsaturatedFlow.hpp.

                                                                             {
    MoFEMFunctionBegin;
    int block_id = -1;
    if (lastEvalBcFluxEnt == ent) {
      block_id = lastEvalBcBlockFluxId;
    } else {
      for (BcMap::iterator it = bcFluxMap.begin(); it != bcFluxMap.end();
           it++) {
        if (it->second->eNts.find(ent) != it->second->eNts.end()) {
          block_id = it->first;
        }
      }
      lastEvalBcFluxEnt = ent;
      lastEvalBcBlockFluxId = block_id;
    }
    if (block_id >= 0) {
      if (bcFluxMap.at(block_id)->hookFun) {
        flux = bcFluxMap.at(block_id)->hookFun(x, y, z);
      } else {
        flux = bcFluxMap.at(block_id)->fixValue;
      }
    } else {
      flux = 0;
    }
    MoFEMFunctionReturn(0);
  }

◆ getBcOnValues()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::getBcOnValues	(	const EntityHandle	ent,
		const int	gg,
		const double	x,
		const double	y,
		const double	z,
		double &	value )

inlinevirtual

Get value on boundary.

Parameters

ent	entity handle
gg	number of integration point
x	x-coordinate
y	y-coordinate
z	z-coordinate
value	returned value

Returns: error code

Reimplemented from MixTransport::MixTransportElement.

Examples: UnsaturatedFlow.hpp.

Definition at line 168 of file UnsaturatedFlow.hpp.

                                              {
    MoFEMFunctionBegin;
    int block_id = -1;
    if (lastEvalBcValEnt == ent) {
      block_id = lastEvalBcBlockValId;
    } else {
      for (BcMap::iterator it = bcValueMap.begin(); it != bcValueMap.end();
           it++) {
        if (it->second->eNts.find(ent) != it->second->eNts.end()) {
          block_id = it->first;
        }
      }
      lastEvalBcValEnt = ent;
      lastEvalBcBlockValId = block_id;
    }
    if (block_id >= 0) {
      if (bcValueMap.at(block_id)->hookFun) {
        value = bcValueMap.at(block_id)->hookFun(x, y, z);
      } else {
        value = bcValueMap.at(block_id)->fixValue;
      }
    } else {
      value = 0;
    }
    MoFEMFunctionReturn(0);
  }

◆ getMaterial()

virtual MoFEMErrorCode MixTransport::UnsaturatedFlowElement::getMaterial	(	const EntityHandle	ent,
		int &	block_id ) const

inlinevirtual

For given element handle get material block Id.

Parameters

ent	finite element entity handle
block_id	reference to returned block id

Returns: error code

Examples: UnsaturatedFlow.hpp.

Definition at line 127 of file UnsaturatedFlow.hpp.

                                                          {
    MoFEMFunctionBegin;
    for (MaterialsDoubleMap::const_iterator mit = dMatMap.begin();
         mit != dMatMap.end(); mit++) {
      if (mit->second->tEts.find(ent) != mit->second->tEts.end()) {
        block_id = mit->first;
        MoFEMFunctionReturnHot(0);
      }
    }
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Element not found, no material data");
    MoFEMFunctionReturn(0);
  }

◆ setFiniteElements()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::setFiniteElements	(	ForcesAndSourcesCore::RuleHookFun	vol_rule = VolRule(),
		ForcesAndSourcesCore::RuleHookFun	face_rule = FaceRule() )

inline

Create finite element instances.

Parameters

vol_rule	integration rule for volume element
face_rule	integration rule for boundary element

Returns: error code

Examples: UnsaturatedFlow.hpp, and unsaturated_transport.cpp.

Definition at line 1446 of file UnsaturatedFlow.hpp.

                                                                            {
    MoFEMFunctionBegin;
 
    // create finite element instances
    feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
        new FaceElementForcesAndSourcesCore(mField));
    feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
        new FaceElementForcesAndSourcesCore(mField));
    feVolInitialPc = boost::shared_ptr<ForcesAndSourcesCore>(
        new VolumeElementForcesAndSourcesCore(mField));
    feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
        new VolumeElementForcesAndSourcesCore(mField));
    feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
        new VolumeElementForcesAndSourcesCore(mField));
    // set integration rule to elements
    feFaceBc->getRuleHook = face_rule;
    feFaceRhs->getRuleHook = face_rule;
    feVolInitialPc->getRuleHook = vol_rule;
    feVolLhs->getRuleHook = vol_rule;
    feVolRhs->getRuleHook = vol_rule;
    // set function hook for finite element postprocessing stage
    feVolRhs->preProcessHook = preProcessVol(*this, feVolRhs);
    feVolLhs->preProcessHook = preProcessVol(*this, feVolLhs);
    feVolRhs->postProcessHook = postProcessVol(*this, feVolRhs);
    feVolLhs->postProcessHook = postProcessVol(*this, feVolLhs);
 
    // Set Piola transform
    CHKERR AddHOOps<2, 3, 3>::add(feFaceBc->getOpPtrVector(), {HDIV});
    CHKERR AddHOOps<2, 3, 3>::add(feFaceRhs->getOpPtrVector(), {HDIV});
 
    // create method for setting history for fluxes on boundary
    scaleMethodFlux = boost::shared_ptr<MethodForForceScaling>(
        new TimeForceScale("-flux_history", false));
 
    // create method for setting history for presure heads on boundary
    scaleMethodValue = boost::shared_ptr<MethodForForceScaling>(
        new TimeForceScale("-head_history", false));
 
    // Set operator to calculate essential boundary conditions
    feFaceBc->getOpPtrVector().push_back(
        new OpEvaluateBcOnFluxes(*this, "FLUXES"));
 
    // Set operator to calculate initial capillary pressure
    CHKERR AddHOOps<3, 3, 3>::add(feVolInitialPc->getOpPtrVector(), {HDIV, L2});
    feVolInitialPc->getOpPtrVector().push_back(
        new OpEvaluateInitiallHead(*this, "VALUES"));
 
    // set residual face from Neumann terms, i.e. applied pressure
    feFaceRhs->getOpPtrVector().push_back(
        new OpRhsBcOnValues(*this, "FLUXES", scaleMethodValue));
    // set residual finite element operators
    headRateAtGaussPts = boost::make_shared<VectorDouble>();
 
 
    // resAtGaussPts = boost::make_shared<VectorDouble>();
    CHKERR AddHOOps<3, 3, 3>::add(feVolRhs->getOpPtrVector(), {HDIV, L2});
    feVolRhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
        string("VALUES") + "_t", headRateAtGaussPts, MBTET));
    feVolRhs->getOpPtrVector().push_back(
        new OpValuesAtGaussPts(*this, "VALUES"));
    feVolRhs->getOpPtrVector().push_back(
        new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
    feVolRhs->getOpPtrVector().push_back(new OpResidualFlux(*this, "FLUXES"));
    feVolRhs->getOpPtrVector().push_back(new OpResidualMass(*this, "VALUES"));
    feVolRhs->getOpPtrVector().back().opType =
        ForcesAndSourcesCore::UserDataOperator::OPROW;
    // set tangent matrix finite element operators
    CHKERR AddHOOps<3, 3, 3>::add(feVolLhs->getOpPtrVector(), {HDIV, L2});
    feVolLhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
        string("VALUES") + "_t", headRateAtGaussPts, MBTET));
    // feVolLhs->getOpPtrVector().push_back(
    //   new
    //   OpCalculateScalarFieldValues(string("FLUXES")+"_residual",resAtGaussPts,MBTET)
    // );
    feVolLhs->getOpPtrVector().push_back(
        new OpValuesAtGaussPts(*this, "VALUES"));
    feVolLhs->getOpPtrVector().push_back(
        new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
    feVolLhs->getOpPtrVector().push_back(
        new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
    feVolLhs->getOpPtrVector().push_back(new OpVU_L2L2(*this, "VALUES"));
    feVolLhs->getOpPtrVector().push_back(
        new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
    feVolLhs->getOpPtrVector().push_back(
        new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
 
    // Adding finite elements to DM, time solver will ask for it to assemble
    // tangent matrices and residuals
    boost::shared_ptr<FEMethod> null;
    CHKERR DMMoFEMTSSetIFunction(dM, "MIX_BCVALUE", feFaceRhs, null, null);
    CHKERR DMMoFEMTSSetIFunction(dM, "MIX", feVolRhs, null, null);
    CHKERR DMMoFEMTSSetIJacobian(dM, "MIX", feVolLhs, null, null);
 
    // setting up post-processing
 
 
    auto get_post_process_ele = [&]() {
      auto post_process = boost::make_shared<PostProcEle>(mField);
 
      CHKERR AddHOOps<3, 3, 3>::add(post_process->getOpPtrVector(), {HDIV, L2});
 
      auto values_ptr = boost::make_shared<VectorDouble>();
      auto grad_ptr = boost::make_shared<MatrixDouble>();
      auto flux_ptr = boost::make_shared<MatrixDouble>();
 
      post_process->getOpPtrVector().push_back(
          new OpCalculateScalarFieldValues("VALUES", values_ptr));
      post_process->getOpPtrVector().push_back(
          new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
      post_process->getOpPtrVector().push_back(
          new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
 
      using OpPPMap = OpPostProcMapInMoab<3, 3>;
 
      post_process->getOpPtrVector().push_back(
 
          new OpPPMap(post_process->getPostProcMesh(),
                      post_process->getMapGaussPts(),
 
                      {{"VALUES", values_ptr}},
 
                      {{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
 
                      {}, {}
 
                      ));
 
      return post_process;
    };
 
    auto post_process = get_post_process_ele();
 
    post_process->getOpPtrVector().push_back(
        new OpValuesAtGaussPts(*this, "VALUES"));
    post_process->getOpPtrVector().push_back(
        new OpPostProcMaterial(*this, post_process->getPostProcMesh(),
                               post_process->getMapGaussPts(), "VALUES"));
 
    // Integrate fluxes on boundary
    boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
    flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
        new FaceElementForcesAndSourcesCore(mField));
    CHKERR AddHOOps<2, 3, 3>::add(flux_integrate->getOpPtrVector(), {HDIV});
    flux_integrate->getOpPtrVector().push_back(
        new OpIntegrateFluxes(*this, "FLUXES"));
    int frequency = 1;
    PetscOptionsBegin(PETSC_COMM_WORLD, "", "Monitor post proc", "none");
    CHKERR PetscOptionsInt(
        "-how_often_output",
        "frequency how often results are dumped on hard disk", "", frequency,
        &frequency, NULL);
    PetscOptionsEnd();
 
    tsMonitor = boost::shared_ptr<FEMethod>(
        new MonitorPostProc(*this, post_process, flux_integrate, frequency));
    TsCtx *ts_ctx;
    DMMoFEMGetTsCtx(dM, &ts_ctx);
    ts_ctx->getPostProcessMonitor().push_back(tsMonitor);
    MoFEMFunctionReturn(0);
  }

◆ solveProblem()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::solveProblem ( bool set_initial_pc = true )

inline

solve problem

Returns: error code

Examples: unsaturated_transport.cpp.

Definition at line 1699 of file UnsaturatedFlow.hpp.

                                                          {
    MoFEMFunctionBegin;
    if (set_initial_pc) {
      // Set initial head
      CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
    }
 
    // Initiate vector from data on the mesh
    CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_FORWARD);
 
    // Create time solver
    TS ts;
    CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
    // Use backward Euler method
    CHKERR TSSetType(ts, TSBEULER);
    // Set final time
    double ftime = 1;
    CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
    // Setup solver from commabd line
    CHKERR TSSetFromOptions(ts);
    // Set DM to TS
    CHKERR TSSetDM(ts, dM);
#if PETSC_VERSION_GE(3, 7, 0)
    CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
#endif
    // Set-up monitor
    TsCtx *ts_ctx;
    DMMoFEMGetTsCtx(dM, &ts_ctx);
    CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULLPTR);
 
    // This add SNES monitor, to show error by fields. It is dirty trick
    // to add monitor, so code is hidden from doxygen
    CHKERR TSSetSolution(ts, D);
    CHKERR TSSetUp(ts);
    SNES snes;
    CHKERR TSGetSNES(ts, &snes);
 
#if PETSC_VERSION_GE(3, 7, 0)
    {
      PetscViewerAndFormat *vf;
      CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
                                        PETSC_VIEWER_DEFAULT, &vf);
      CHKERR SNESMonitorSet(
          snes,
          (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
                             void *))SNESMonitorFields,
          vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
    }
#else
    {
      CHKERR SNESMonitorSet(snes,
                            (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
                                               void *))SNESMonitorFields,
                            0, 0);
    }
#endif
 
    CHKERR TSSolve(ts, D);
 
    // Get statistic form TS and print it
    CHKERR TSGetTime(ts, &ftime);
    PetscInt steps, snesfails, rejects, nonlinits, linits;
    CHKERR TSGetTimeStepNumber(ts, &steps);
    CHKERR TSGetSNESFailures(ts, &snesfails);
    CHKERR TSGetStepRejections(ts, &rejects);
    CHKERR TSGetSNESIterations(ts, &nonlinits);
    CHKERR TSGetKSPIterations(ts, &linits);
    PetscPrintf(PETSC_COMM_WORLD,
                "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
                "%D, linits %D\n",
                steps, rejects, snesfails, ftime, nonlinits, linits);
 
    MoFEMFunctionReturn(0);
  }