Mix transport problem. More...

#include <users_modules/tutorials/cor-0to1/src/MixTransportElement.hpp>

Inheritance diagram for MixTransport::MixTransportElement:

Collaboration diagram for MixTransport::MixTransportElement:

Classes
struct	BlockData
	data for evaluation of het conductivity and heat capacity elements More...

struct	MyTriFE
	definition of surface element More...

struct	MyVolumeFE
	definition of volume element More...

struct	OpDivTauU_HdivL2
	Assemble \( \int_\mathcal{T} u \textrm{div}[\boldsymbol\tau] \textrm{d}\mathcal{T} \). More...

struct	OpError
	calculate error evaluator More...

struct	OpEvaluateBcOnFluxes
	Evaluate boundary conditions on fluxes. More...

struct	OpFluxDivergenceAtGaussPts
	calculate flux at integration point More...

struct	OpL2Source
	Calculate source therms, i.e. \(\int_\mathcal{T} f v \textrm{d}\mathcal{T}\). More...

struct	OpPostProc

struct	OpRhsBcOnValues
	calculate \( \int_\mathcal{S} {\boldsymbol\tau} \cdot \mathbf{n}u \textrm{d}\mathcal{S} \) More...

struct	OpSkeleton
	calculate jump on entities More...

struct	OpTauDotSigma_HdivHdiv
	Assemble \(\int_\mathcal{T} \mathbf{A} \boldsymbol\sigma \cdot \boldsymbol\tau \textrm{d}\mathcal{T}\). More...

struct	OpValuesAtGaussPts
	Calculate values at integration points. More...

struct	OpValuesGradientAtGaussPts
	Calculate gradients of values at integration points. More...

struct	OpVDivSigma_L2Hdiv
	\( \int_\mathcal{T} \textrm{div}[\boldsymbol\sigma] v \textrm{d}\mathcal{T} \) More...

Public Member Functions
	MixTransportElement (MoFEM::Interface &m_field)
	construction of this data structure More...

virtual	~MixTransportElement ()
	destructor More...

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

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

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

virtual MoFEMErrorCode	getBcOnValues (const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
	evaluate natural (Dirichlet) boundary conditions More...

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

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

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

MoFEMErrorCode	buildProblem (BitRefLevel &ref_level)
	Build problem. More...

MoFEMErrorCode	postProc (const string out_file)
	Post process results. More...

MoFEMErrorCode	createMatrices ()
	create matrices More...

MoFEMErrorCode	solveLinearProblem ()
	solve problem More...

MoFEMErrorCode	calculateResidual ()
	calculate residual More...

MoFEMErrorCode	evaluateError ()
	Calculate error on elements. More...

MoFEMErrorCode	destroyMatrices ()
	destroy matrices More...

Public Attributes
MoFEM::Interface &	mField

MyVolumeFE	feVol
	Instance of volume element. More...

MyTriFE	feTri
	Instance of surface element. More...

VectorDouble	valuesAtGaussPts
	values at integration points on element More...

MatrixDouble	valuesGradientAtGaussPts
	gradients at integration points on element More...

VectorDouble	divergenceAtGaussPts
	divergence at integration points on element More...

MatrixDouble	fluxesAtGaussPts
	fluxes at integration points on element More...

set< int >	bcIndices

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

Vec	D

Vec	D0

Vec	F

Mat	Aij

map< double, EntityHandle >	errorMap

double	sumErrorFlux

double	sumErrorDiv

double	sumErrorJump

Detailed Description

Mix transport problem.

Note to solve this system you need to use direct solver or proper preconditioner for saddle problem.

It is based on [6] [7] [reviartthomas1996] https://www.researchgate.net/profile/Richard_Falk/publication/226454406_Differential_Complexes_and_Stability_of_Finite_Element_Methods_I._The_de_Rham_Complex/links/02e7e5214f0426ff77000000.pdf

General problem have form,

\[ \mathbf{A} \boldsymbol\sigma + \textrm{grad}[u] = \mathbf{0} \; \textrm{on} \; \Omega \\ \textrm{div}[\boldsymbol\sigma] = f \; \textrm{on} \; \Omega \]

Definition at line 31 of file MixTransportElement.hpp.

Constructor & Destructor Documentation

◆ MixTransportElement()

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

construction of this data structure

Examples: UnsaturatedFlow.hpp.

Definition at line 67 of file MixTransportElement.hpp.

68 : mField(m_field), feVol(m_field), feTri(m_field){};

MixTransport::MixTransportElement::mField

MoFEM::Interface & mField

Definition: MixTransportElement.hpp:33

MixTransport::MixTransportElement::feTri

MyTriFE feTri

Instance of surface element.

Definition: MixTransportElement.hpp:62

MixTransport::MixTransportElement::feVol

MyVolumeFE feVol

Instance of volume element.

Definition: MixTransportElement.hpp:48

◆ ~MixTransportElement()

virtual MixTransport::MixTransportElement::~MixTransportElement ( )

virtual

destructor

Definition at line 73 of file MixTransportElement.hpp.

73{}

Member Function Documentation

◆ addFields()

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

Add fields to database.

Parameters

values	name of the fields
fluxes	name of filed for fluxes
order	order of approximation

Returns: error code

Definition at line 191 of file MixTransportElement.hpp.

                                            {
    MoFEMFunctionBegin;
    // Fields
    CHKERR mField.add_field(fluxes, HDIV, DEMKOWICZ_JACOBI_BASE, 1);
    CHKERR mField.add_field(values, L2, AINSWORTH_LEGENDRE_BASE, 1);
 
    // meshset consisting all entities in mesh
    EntityHandle root_set = mField.get_moab().get_root_set();
    // add entities to field
    CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, fluxes);
    CHKERR mField.add_ents_to_field_by_type(root_set, MBTET, values);
    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);
    MoFEMFunctionReturn(0);
  }

◆ addFiniteElements()

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

add finite elements

Definition at line 210 of file MixTransportElement.hpp.

                                                                  {
    MoFEMFunctionBegin;
 
    // Set up volume element operators. Operators are used to calculate
    // components of stiffness matrix & right hand side, in essence are used to
    // do volume integrals over tetrahedral in this case.
 
    // 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);
 
    // Define finite element to integrate over skeleton, we need that to
    // evaluate error
    CHKERR mField.add_finite_element("MIX_SKELETON", MF_ZERO);
    CHKERR mField.modify_finite_element_add_field_row("MIX_SKELETON",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_col("MIX_SKELETON",
                                                      fluxes_name);
    CHKERR mField.modify_finite_element_add_field_data("MIX_SKELETON",
                                                       fluxes_name);
 
    // Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those
    // BLOCKSETS and add them to "MIX" finite element. In addition get data form
    // that meshset and set conductivity which is used to calculate fluxes from
    // gradients of concentration or gradient of temperature, depending how you
    // interpret variables.
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             mField, BLOCKSET | MAT_THERMALSET, it)) {
 
      Mat_Thermal temp_data;
      CHKERR it->getAttributeDataStructure(temp_data);
      setOfBlocks[it->getMeshsetId()].cOnductivity =
          temp_data.data.Conductivity;
      setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
      CHKERR mField.get_moab().get_entities_by_type(
          it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
      CHKERR mField.add_ents_to_finite_element_by_type(
          setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
 
      Range skeleton;
      CHKERR mField.get_moab().get_adjacencies(
          setOfBlocks[it->getMeshsetId()].tEts, 2, false, skeleton,
          moab::Interface::UNION);
      CHKERR mField.add_ents_to_finite_element_by_type(skeleton, MBTRI,
                                                       "MIX_SKELETON");
    }
 
    // 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);
 
    // 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);
    MoFEMFunctionReturn(0);
  }

◆ buildProblem()

MoFEMErrorCode MixTransport::MixTransportElement::buildProblem ( BitRefLevel & ref_level )

Build problem.

Parameters

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

Returns: error code

Definition at line 293 of file MixTransportElement.hpp.

                                                      {
    MoFEMFunctionBegin;
    // build field
    CHKERR mField.build_fields();
    // get tetrahedrons which has been build previously and now in so called
    // garbage bit level
    Range done_tets;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        BitRefLevel().set(0), BitRefLevel().set(), MBTET, done_tets);
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTET,
        done_tets);
    // get tetrahedrons which belong to problem bit level
    Range ref_tets;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        ref_level, BitRefLevel().set(), MBTET, ref_tets);
    ref_tets = subtract(ref_tets, done_tets);
    CHKERR mField.build_finite_elements("MIX", &ref_tets, 2);
    // get triangles which has been build previously and now in so called
    // garbage bit level
    Range done_faces;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        BitRefLevel().set(0), BitRefLevel().set(), MBTRI, done_faces);
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTRI,
        done_faces);
    // get triangles which belong to problem bit level
    Range ref_faces;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        ref_level, BitRefLevel().set(), MBTRI, ref_faces);
    ref_faces = subtract(ref_faces, done_faces);
    // build finite elements structures
    CHKERR mField.build_finite_elements("MIX_BCFLUX", &ref_faces, 2);
    CHKERR mField.build_finite_elements("MIX_BCVALUE", &ref_faces, 2);
    CHKERR mField.build_finite_elements("MIX_SKELETON", &ref_faces, 2);
    // Build adjacencies of degrees of freedom and elements
    CHKERR mField.build_adjacencies(ref_level);
    // Define problem
    CHKERR mField.add_problem("MIX", MF_ZERO);
    // set refinement level for problem
    CHKERR mField.modify_problem_ref_level_set_bit("MIX", ref_level);
    // Add element to problem
    CHKERR mField.modify_problem_add_finite_element("MIX", "MIX");
    CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_SKELETON");
    // Boundary conditions
    CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCFLUX");
    CHKERR mField.modify_problem_add_finite_element("MIX", "MIX_BCVALUE");
    // build problem
 
    ProblemsManager *prb_mng_ptr;
    CHKERR mField.getInterface(prb_mng_ptr);
    CHKERR prb_mng_ptr->buildProblem("MIX", true);
    // mesh partitioning
    // partition
    CHKERR prb_mng_ptr->partitionProblem("MIX");
    CHKERR prb_mng_ptr->partitionFiniteElements("MIX");
    // what are ghost nodes, see Petsc Manual
    CHKERR prb_mng_ptr->partitionGhostDofs("MIX");
    MoFEMFunctionReturn(0);
  }

◆ calculateResidual()

MoFEMErrorCode MixTransport::MixTransportElement::calculateResidual ( )

calculate residual

Definition at line 624 of file MixTransportElement.hpp.

                                     {
    MoFEMFunctionBegin;
    CHKERR VecZeroEntries(F);
    CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecAssemblyBegin(F);
    CHKERR VecAssemblyEnd(F);
    // calculate residuals
    feVol.getOpPtrVector().clear();
    feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
    feVol.getOpPtrVector().push_back(
        new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
    feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
    feVol.getOpPtrVector().push_back(
        new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
    feVol.getOpPtrVector().push_back(
        new OpTauDotSigma_HdivHdiv(*this, "FLUXES", PETSC_NULL, F));
    feVol.getOpPtrVector().push_back(
        new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", PETSC_NULL, F));
    CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
    feTri.getOpPtrVector().clear();
    feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
    CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
    CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(F);
    CHKERR VecAssemblyEnd(F);
    // CHKERR VecAXPY(F,-1.,D0);
    // CHKERR MatZeroRowsIS(Aij,essential_bc_ids,1,PETSC_NULL,F);
    {
      std::vector<int> ids;
      ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
      std::vector<double> vals(ids.size(), 0);
      CHKERR VecSetValues(F, ids.size(), &*ids.begin(), &*vals.begin(),
                          INSERT_VALUES);
      CHKERR VecAssemblyBegin(F);
      CHKERR VecAssemblyEnd(F);
    }
    {
      double nrm2_F;
      CHKERR VecNorm(F, NORM_2, &nrm2_F);
      PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
      const double eps = 1e-8;
      if (nrm2_F > eps) {
        // SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"problem with
        // residual");
      }
    }
    MoFEMFunctionReturn(0);
  }

◆ createMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::createMatrices ( )

create matrices

Definition at line 485 of file MixTransportElement.hpp.

                                  {
    MoFEMFunctionBegin;
    CHKERR mField.getInterface<MatrixManager>()
        ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("MIX", &Aij);
    CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D);
    CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", COL, &D0);
    CHKERR mField.getInterface<VecManager>()->vecCreateGhost("MIX", ROW, &F);
    MoFEMFunctionReturn(0);
  }

◆ destroyMatrices()

MoFEMErrorCode MixTransport::MixTransportElement::destroyMatrices ( )

destroy matrices

Definition at line 710 of file MixTransportElement.hpp.

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

◆ evaluateError()

MoFEMErrorCode MixTransport::MixTransportElement::evaluateError ( )

Calculate error on elements.

Todo:: this functions runs serial, in future should be parallel and work on distributed meshes.

Definition at line 681 of file MixTransportElement.hpp.

                                 {
    MoFEMFunctionBegin;
    errorMap.clear();
    sumErrorFlux = 0;
    sumErrorDiv = 0;
    sumErrorJump = 0;
    feTri.getOpPtrVector().clear();
    feTri.getOpPtrVector().push_back(new OpSkeleton(*this, "FLUXES"));
    CHKERR mField.loop_finite_elements("MIX", "MIX_SKELETON", feTri, 0,
                                       mField.get_comm_size());
    feVol.getOpPtrVector().clear();
    feVol.getOpPtrVector().push_back(
        new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
    feVol.getOpPtrVector().push_back(
        new OpValuesGradientAtGaussPts(*this, "VALUES"));
    feVol.getOpPtrVector().push_back(new OpError(*this, "VALUES"));
    CHKERR mField.loop_finite_elements("MIX", "MIX", feVol, 0,
                                       mField.get_comm_size());
    const Problem *problem_ptr;
    CHKERR mField.get_problem("MIX", &problem_ptr);
    PetscPrintf(mField.get_comm(),
                "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
                "jump^2 = %6.4e error tot^2 = %6.4e\n",
                problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
                sumErrorJump, sumErrorFlux + sumErrorDiv + sumErrorJump);
    MoFEMFunctionReturn(0);
  }

◆ getBcOnFluxes()

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

virtual

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 in MyTransport, MyTransport, and MixTransport::UnsaturatedFlowElement.

Definition at line 165 of file MixTransportElement.hpp.

                                                     {
    MoFEMFunctionBeginHot;
    flux = 0;
    MoFEMFunctionReturnHot(0);
  }

◆ getBcOnValues()

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

virtual

evaluate natural (Dirichlet) boundary conditions

Parameters

ent	entity on which bc is evaluated
x	coordinate
y	coordinate
z	coordinate
value	vale

Returns: error code

Reimplemented in MixTransport::UnsaturatedFlowElement.

Definition at line 148 of file MixTransportElement.hpp.

                                                                      {
    MoFEMFunctionBeginHot;
    value = 0;
    MoFEMFunctionReturnHot(0);
  }

◆ getDirichletBCIndices()

MoFEMErrorCode MixTransport::MixTransportElement::getDirichletBCIndices ( IS * is )

get dof indices where essential boundary conditions are applied

Parameters

is indices

Returns: error code

Definition at line 89 of file MixTransportElement.hpp.

                                               {
    MoFEMFunctionBegin;
    std::vector<int> ids;
    ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
    IS is_local;
    CHKERR ISCreateGeneral(mField.get_comm(), ids.size(),
                           ids.empty() ? PETSC_NULL : &ids[0],
                           PETSC_COPY_VALUES, &is_local);
    CHKERR ISAllGather(is_local, is);
    CHKERR ISDestroy(&is_local);
    MoFEMFunctionReturn(0);
  }

◆ getResistivity()

virtual MoFEMErrorCode MixTransport::MixTransportElement::getResistivity	(	const EntityHandle	ent,
		const double	x,
		const double	y,
		const double	z,
		MatrixDouble3by3 &	inv_k
	)

virtual

natural (Dirichlet) boundary conditions (set values)

Parameters

ent	handle to finite element entity
x	coord
y	coord
z	coord
value	reference to value set by function

Returns: error code

Definition at line 128 of file MixTransportElement.hpp.

                                                                 {
    MoFEMFunctionBeginHot;
    inv_k.clear();
    for (int dd = 0; dd < 3; dd++) {
      inv_k(dd, dd) = 1;
    }
    MoFEMFunctionReturnHot(0);
  }

◆ getSource()

virtual MoFEMErrorCode MixTransport::MixTransportElement::getSource	(	const EntityHandle	ent,
		const double	x,
		const double	y,
		const double	z,
		double &	flux
	)

virtual

set source term

Parameters

ent	handle to entity on which function is evaluated
x	coord
y	coord
z	coord
flux	reference to source term set by function

Returns: error code

Reimplemented in MyTransport, and MyTransport.

Definition at line 111 of file MixTransportElement.hpp.

                                                 {
    MoFEMFunctionBeginHot;
    flux = 0;
    MoFEMFunctionReturnHot(0);
  }

◆ postProc()

MoFEMErrorCode MixTransport::MixTransportElement::postProc ( const string out_file )

Post process results.

Returns: error code

Definition at line 430 of file MixTransportElement.hpp.

                                                 {
    MoFEMFunctionBegin;
 
    PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
        mField);
 
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    post_proc.getOpPtrVector().push_back(new OpCalculateHOJac<3>(jac_ptr));
    post_proc.getOpPtrVector().push_back(
        new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
    post_proc.getOpPtrVector().push_back(
        new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
    post_proc.getOpPtrVector().push_back(
        new OpSetHOInvJacToScalarBases<3>(L2, inv_jac_ptr));
 
    auto values_ptr = boost::make_shared<VectorDouble>();
    auto grad_ptr = boost::make_shared<MatrixDouble>();
    auto flux_ptr = boost::make_shared<MatrixDouble>();
 
    post_proc.getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("VALUES", values_ptr));
    post_proc.getOpPtrVector().push_back(
        new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
    post_proc.getOpPtrVector().push_back(
        new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
 
    using OpPPMap = OpPostProcMapInMoab<3, 3>;
 
    post_proc.getOpPtrVector().push_back(
 
        new OpPPMap(post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
 
                    {{"VALUES", values_ptr}},
 
                    {{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
 
                    {}, {}
 
                    ));
 
    post_proc.getOpPtrVector().push_back(new OpPostProc(
        post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
 
    CHKERR mField.loop_finite_elements("MIX", "MIX", post_proc);
 
    CHKERR post_proc.writeFile(out_file.c_str());
    MoFEMFunctionReturn(0);
  }

◆ solveLinearProblem()

MoFEMErrorCode MixTransport::MixTransportElement::solveLinearProblem ( )

solve problem

Returns: error code

Definition at line 499 of file MixTransportElement.hpp.

                                      {
    MoFEMFunctionBegin;
    CHKERR MatZeroEntries(Aij);
    CHKERR VecZeroEntries(F);
    CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecZeroEntries(D0);
    CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecZeroEntries(D);
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
 
    CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
        "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
 
    // Calculate essential boundary conditions
 
    // clear essential bc indices, it could have dofs from other mesh refinement
    bcIndices.clear();
    // clear operator, just in case if some other operators are left on this
    // element
    feTri.getOpPtrVector().clear();
    feTri.getOpPtrVector().push_back(
        new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
    // set operator to calculate essential boundary conditions
    feTri.getOpPtrVector().push_back(new OpEvaluateBcOnFluxes(*this, "FLUXES"));
    CHKERR mField.loop_finite_elements("MIX", "MIX_BCFLUX", feTri);
    CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(D0);
    CHKERR VecAssemblyEnd(D0);
 
    // set operators to calculate matrix and right hand side vectors
    feVol.getOpPtrVector().clear();
    feVol.getOpPtrVector().push_back(new OpL2Source(*this, "VALUES", F));
    feVol.getOpPtrVector().push_back(
        new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
    feVol.getOpPtrVector().push_back(new OpValuesAtGaussPts(*this, "VALUES"));
    feVol.getOpPtrVector().push_back(
        new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES", F));
    feVol.getOpPtrVector().push_back(
        new OpTauDotSigma_HdivHdiv(*this, "FLUXES", Aij, F));
    feVol.getOpPtrVector().push_back(
        new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES", Aij, F));
    CHKERR mField.loop_finite_elements("MIX", "MIX", feVol);
 
    // calculate right hand side for natural boundary conditions
    feTri.getOpPtrVector().clear();
    feTri.getOpPtrVector().push_back(
        new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
    feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this, "FLUXES", F));
    CHKERR mField.loop_finite_elements("MIX", "MIX_BCVALUE", feTri);
 
    // assemble matrices
    CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
    CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(F);
    CHKERR VecAssemblyEnd(F);
 
    {
      double nrm2_F;
      CHKERR VecNorm(F, NORM_2, &nrm2_F);
      PetscPrintf(PETSC_COMM_WORLD, "nrm2_F = %6.4e\n", nrm2_F);
    }
 
    // CHKERR MatMultAdd(Aij,D0,F,F); ;
 
    // for ksp solver vector is moved into rhs side
    // for snes it is left ond the left
    CHKERR VecScale(F, -1);
 
    IS essential_bc_ids;
    CHKERR getDirichletBCIndices(&essential_bc_ids);
    CHKERR MatZeroRowsColumnsIS(Aij, essential_bc_ids, 1, D0, F);
    CHKERR ISDestroy(&essential_bc_ids);
 
    // {
    //   double norm;
    //   CHKERR MatNorm(Aij,NORM_FROBENIUS,&norm); ;
    //   PetscPrintf(PETSC_COMM_WORLD,"mat norm = %6.4e\n",norm);
    // }
 
    {
      double nrm2_F;
      CHKERR VecNorm(F, NORM_2, &nrm2_F);
      PetscPrintf(PETSC_COMM_WORLD, "With BC nrm2_F = %6.4e\n", nrm2_F);
    }
 
    // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
    // MatView(Aij,PETSC_VIEWER_STDOUT_WORLD);
    // std::string wait;
    // std::cin >> wait;
 
    // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);
    // std::cin >> wait;
 
    // Solve
    KSP solver;
    CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
    CHKERR KSPSetOperators(solver, Aij, Aij);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetUp(solver);
    CHKERR KSPSolve(solver, F, D);
    CHKERR KSPDestroy(&solver);
 
    {
      double nrm2_D;
      CHKERR VecNorm(D, NORM_2, &nrm2_D);
      ;
      PetscPrintf(PETSC_COMM_WORLD, "nrm2_D = %6.4e\n", nrm2_D);
    }
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
 
    // copy data form vector on mesh
    CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
        "MIX", COL, D, INSERT_VALUES, SCATTER_REVERSE);
 
    MoFEMFunctionReturn(0);
  }

Member Data Documentation

◆ Aij

Mat MixTransport::MixTransportElement::Aij

Definition at line 482 of file MixTransportElement.hpp.

◆ bcIndices

set<int> MixTransport::MixTransportElement::bcIndices

Examples: UnsaturatedFlow.hpp.

Definition at line 82 of file MixTransportElement.hpp.

◆ D

Vec MixTransport::MixTransportElement::D

Definition at line 481 of file MixTransportElement.hpp.

◆ D0

Vec MixTransport::MixTransportElement::D0

Examples: UnsaturatedFlow.hpp.

Definition at line 481 of file MixTransportElement.hpp.

◆ divergenceAtGaussPts

VectorDouble MixTransport::MixTransportElement::divergenceAtGaussPts

divergence at integration points on element

Examples: UnsaturatedFlow.hpp.

Definition at line 79 of file MixTransportElement.hpp.

◆ errorMap

map<double, EntityHandle> MixTransport::MixTransportElement::errorMap

Definition at line 1370 of file MixTransportElement.hpp.

◆ F

Vec MixTransport::MixTransportElement::F

Definition at line 481 of file MixTransportElement.hpp.

◆ feTri

MyTriFE MixTransport::MixTransportElement::feTri

Instance of surface element.

Definition at line 62 of file MixTransportElement.hpp.

◆ feVol

MyVolumeFE MixTransport::MixTransportElement::feVol

Instance of volume element.

Definition at line 48 of file MixTransportElement.hpp.

◆ fluxesAtGaussPts

MatrixDouble MixTransport::MixTransportElement::fluxesAtGaussPts

fluxes at integration points on element

Examples: UnsaturatedFlow.hpp.

Definition at line 80 of file MixTransportElement.hpp.

◆ mField

MoFEM::Interface& MixTransport::MixTransportElement::mField

Examples: UnsaturatedFlow.hpp.

Definition at line 33 of file MixTransportElement.hpp.

◆ setOfBlocks

std::map<int, BlockData> MixTransport::MixTransportElement::setOfBlocks

maps block set id with appropriate BlockData

Definition at line 182 of file MixTransportElement.hpp.

◆ sumErrorDiv

double MixTransport::MixTransportElement::sumErrorDiv

Definition at line 1372 of file MixTransportElement.hpp.

◆ sumErrorFlux

double MixTransport::MixTransportElement::sumErrorFlux

Definition at line 1371 of file MixTransportElement.hpp.

◆ sumErrorJump

double MixTransport::MixTransportElement::sumErrorJump

Definition at line 1373 of file MixTransportElement.hpp.

◆ valuesAtGaussPts

VectorDouble MixTransport::MixTransportElement::valuesAtGaussPts

values at integration points on element

Examples: UnsaturatedFlow.hpp.

Definition at line 75 of file MixTransportElement.hpp.

◆ valuesGradientAtGaussPts

MatrixDouble MixTransport::MixTransportElement::valuesGradientAtGaussPts

gradients at integration points on element

Definition at line 77 of file MixTransportElement.hpp.

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

users_modules/tutorials/cor-0to1/src/MixTransportElement.hpp

Classes

Public Member Functions

Public Attributes

Detailed Description

Constructor & Destructor Documentation

◆ MixTransportElement()

◆ ~MixTransportElement()

Member Function Documentation

◆ addFields()

◆ addFiniteElements()

◆ buildProblem()

◆ calculateResidual()

◆ createMatrices()

◆ destroyMatrices()

◆ evaluateError()

◆ getBcOnFluxes()

◆ getBcOnValues()

◆ getDirichletBCIndices()

◆ getResistivity()

◆ getSource()

◆ postProc()

◆ solveLinearProblem()

Member Data Documentation

◆ Aij

◆ bcIndices

◆ D

◆ D0

◆ divergenceAtGaussPts

◆ errorMap

◆ F

◆ feTri

◆ feVol

◆ fluxesAtGaussPts

◆ mField

◆ setOfBlocks

◆ sumErrorDiv

◆ sumErrorFlux

◆ sumErrorJump

◆ valuesAtGaussPts

◆ valuesGradientAtGaussPts