NonlinearPoisson Struct Reference

Collaboration diagram for NonlinearPoisson:

Public Member Functions
	NonlinearPoisson (MoFEM::Interface &m_field)
MoFEMErrorCode	runProgram ()

Private Types
using	PostProcEle = PostProcBrokenMeshInMoab< FaceEle >

Private Member Functions
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	setIntegrationRules ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()

Static Private Member Functions
static double	sourceTermFunction (const double x, const double y, const double z)
static double	boundaryFunction (const double x, const double y, const double z)

Private Attributes
MoFEM::Interface &	mField
Simple *	simpleInterface
MPI_Comm	mpiComm
const int	mpiRank
SmartPetscObj< DM >	dM
SmartPetscObj< SNES >	snesSolver
std::string	domainField
int	order
boost::shared_ptr< std::vector< unsigned char > >	boundaryMarker
boost::shared_ptr< FaceEle >	domainTangentMatrixPipeline
boost::shared_ptr< FaceEle >	domainResidualVectorPipeline
boost::shared_ptr< EdgeEle >	boundaryTangentMatrixPipeline
boost::shared_ptr< EdgeEle >	boundaryResidualVectorPipeline
boost::shared_ptr< DataAtGaussPoints >	previousUpdate
boost::shared_ptr< VectorDouble >	fieldValuePtr
boost::shared_ptr< MatrixDouble >	fieldGradPtr
boost::shared_ptr< PostProcEle >	postProc
Range	boundaryEntitiesForFieldsplit

Detailed Description

Definition at line 10 of file nonlinear_poisson_2d.cpp.

Member Typedef Documentation

PostProcEle

using NonlinearPoisson::PostProcEle = PostProcBrokenMeshInMoab<FaceEle>

private

Definition at line 72 of file nonlinear_poisson_2d.cpp.

Constructor & Destructor Documentation

NonlinearPoisson()

NonlinearPoisson::NonlinearPoisson

(

MoFEM::Interface &

m_field

)

Definition at line 81 of file nonlinear_poisson_2d.cpp.

    : domainField("U"), mField(m_field), mpiComm(mField.get_comm()),
      mpiRank(mField.get_comm_rank()) {
  domainTangentMatrixPipeline = boost::shared_ptr<FaceEle>(new FaceEle(mField));
  domainResidualVectorPipeline =
      boost::shared_ptr<FaceEle>(new FaceEle(mField));
  boundaryTangentMatrixPipeline =
      boost::shared_ptr<EdgeEle>(new EdgeEle(mField));
  boundaryResidualVectorPipeline =
      boost::shared_ptr<EdgeEle>(new EdgeEle(mField));
 
  previousUpdate =
      boost::shared_ptr<DataAtGaussPoints>(new DataAtGaussPoints());
  fieldValuePtr = boost::shared_ptr<VectorDouble>(previousUpdate,
                                                  &previousUpdate->fieldValue);
  fieldGradPtr = boost::shared_ptr<MatrixDouble>(previousUpdate,
                                                 &previousUpdate->fieldGrad);
}

Member Function Documentation

assembleSystem()

MoFEMErrorCode NonlinearPoisson::assembleSystem ( )

private

Definition at line 198 of file nonlinear_poisson_2d.cpp.

                                                {
  MoFEMFunctionBegin;
 
  auto det_ptr = boost::make_shared<VectorDouble>();
  auto jac_ptr = boost::make_shared<MatrixDouble>();
  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
 
  { // Push operators to the Pipeline that is responsible for calculating the
    // domain tangent matrix (LHS)
 
    // Add default operators to calculate inverse of Jacobian (needed for
    // implementation of 2D problem but not 3D ones)
 
    domainTangentMatrixPipeline->getOpPtrVector().push_back(
        new OpCalculateHOJac<2>(jac_ptr));
    domainTangentMatrixPipeline->getOpPtrVector().push_back(
        new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
    domainTangentMatrixPipeline->getOpPtrVector().push_back(
        new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
    domainTangentMatrixPipeline->getOpPtrVector().push_back(
        new OpSetHOWeightsOnFace());
 
    // Add default operator to calculate field values at integration points
    domainTangentMatrixPipeline->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues(domainField, fieldValuePtr));
    // Add default operator to calculate field gradient at integration points
    domainTangentMatrixPipeline->getOpPtrVector().push_back(
        new OpCalculateScalarFieldGradient<2>(domainField, fieldGradPtr));
 
    // Push operators for domain tangent matrix (LHS)
    domainTangentMatrixPipeline->getOpPtrVector().push_back(
        new OpDomainTangentMatrix(domainField, domainField, previousUpdate,
                                  boundaryMarker));
  }
 
  { // Push operators to the Pipeline that is responsible for calculating the
    // domain residual vector (RHS)
 
    // Add default operators to calculate inverse of Jacobian (needed for
    // implementation of 2D problem but not 3D ones)
    domainResidualVectorPipeline->getOpPtrVector().push_back(
        new OpCalculateHOJac<2>(jac_ptr));
    domainResidualVectorPipeline->getOpPtrVector().push_back(
        new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
    domainResidualVectorPipeline->getOpPtrVector().push_back(
        new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
    domainResidualVectorPipeline->getOpPtrVector().push_back(
        new OpSetHOWeightsOnFace());
 
    // Add default operator to calculate field values at integration points
    domainResidualVectorPipeline->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues(domainField, fieldValuePtr));
    // Add default operator to calculate field gradient at integration points
    domainResidualVectorPipeline->getOpPtrVector().push_back(
        new OpCalculateScalarFieldGradient<2>(domainField, fieldGradPtr));
 
    domainResidualVectorPipeline->getOpPtrVector().push_back(
        new OpDomainResidualVector(domainField, sourceTermFunction,
                                   previousUpdate, boundaryMarker));
  }
 
  { // Push operators to the Pipeline that is responsible for calculating the
    // boundary tangent matrix (LHS)
 
    boundaryTangentMatrixPipeline->getOpPtrVector().push_back(
        new OpBoundaryTangentMatrix(domainField, domainField));
  }
 
  { // Push operators to the Pipeline that is responsible for calculating
    // boundary residual vector (RHS)
 
    // Add default operator to calculate field values at integration points
    boundaryResidualVectorPipeline->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues(domainField, fieldValuePtr));
 
    boundaryResidualVectorPipeline->getOpPtrVector().push_back(
        new OpBoundaryResidualVector(domainField, boundaryFunction,
                                     previousUpdate));
  }
 
  // get Discrete Manager (SmartPetscObj)
  dM = simpleInterface->getDM();
 
  { // Set operators for nonlinear equations solver (SNES) from MoFEM Pipelines
 
    // Set operators for calculation of LHS and RHS of domain elements
    boost::shared_ptr<FaceEle> null_face;
    CHKERR DMMoFEMSNESSetJacobian(dM, simpleInterface->getDomainFEName(),
                                  domainTangentMatrixPipeline, null_face,
                                  null_face);
    CHKERR DMMoFEMSNESSetFunction(dM, simpleInterface->getDomainFEName(),
                                  domainResidualVectorPipeline, null_face,
                                  null_face);
 
    // Set operators for calculation of LHS and RHS of boundary elements
    boost::shared_ptr<EdgeEle> null_edge;
    CHKERR DMMoFEMSNESSetJacobian(dM, simpleInterface->getBoundaryFEName(),
                                  boundaryTangentMatrixPipeline, null_edge,
                                  null_edge);
    CHKERR DMMoFEMSNESSetFunction(dM, simpleInterface->getBoundaryFEName(),
                                  boundaryResidualVectorPipeline, null_edge,
                                  null_edge);
  }
 
  MoFEMFunctionReturn(0);
}

boundaryCondition()

MoFEMErrorCode NonlinearPoisson::boundaryCondition ( )

private

Definition at line 157 of file nonlinear_poisson_2d.cpp.

                                                   {
  MoFEMFunctionBegin;
 
  // Get boundary edges marked in block named "BOUNDARY_CONDITION"
  auto get_ents_on_mesh_skin = [&]() {
    Range boundary_entities;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      std::string entity_name = it->getName();
      if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
        CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
                                                   boundary_entities, true);
      }
    }
    // Add vertices to boundary entities
    Range boundary_vertices;
    CHKERR mField.get_moab().get_connectivity(boundary_entities,
                                              boundary_vertices, true);
    boundary_entities.merge(boundary_vertices);
 
    // Store entities for fieldsplit (block) solver
    boundaryEntitiesForFieldsplit = boundary_entities;
 
    return boundary_entities;
  };
 
  auto mark_boundary_dofs = [&](Range &&skin_edges) {
    auto problem_manager = mField.getInterface<ProblemsManager>();
    auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
    problem_manager->markDofs(simpleInterface->getProblemName(), ROW,
                              skin_edges, *marker_ptr);
    return marker_ptr;
  };
 
  // Get global local vector of marked DOFs. Is global, since is set for all
  // DOFs on processor. Is local since only DOFs on processor are in the
  // vector. To access DOFs use local indices.
  boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
 
  MoFEMFunctionReturn(0);
}

boundaryFunction()

static double NonlinearPoisson::boundaryFunction ( const double  x, const double  y, const double  z  )

inlinestaticprivate

Definition at line 35 of file nonlinear_poisson_2d.cpp.

                                                 {
    return sin(x * 10.) * cos(y * 10.);
    // return 0;
  }

outputResults()

MoFEMErrorCode NonlinearPoisson::outputResults ( )

private

Definition at line 357 of file nonlinear_poisson_2d.cpp.

                                               {
  MoFEMFunctionBegin;
 
  postProc = boost::make_shared<PostProcEle>(mField);
 
  auto u_ptr = boost::make_shared<VectorDouble>();
  postProc->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues(domainField, u_ptr));
 
  using OpPPMap = OpPostProcMapInMoab<2, 2>;
 
  postProc->getOpPtrVector().push_back(
 
      new OpPPMap(postProc->getPostProcMesh(), postProc->getMapGaussPts(),
 
                  {{"U", u_ptr}},
 
                  {}, {}, {}
 
                  )
 
  );
 
  CHKERR DMoFEMLoopFiniteElements(
      dM, simpleInterface->getDomainFEName(),
      boost::dynamic_pointer_cast<FEMethod>(postProc));
 
  CHKERR postProc->writeFile("out_result.h5m");
 
  MoFEMFunctionReturn(0);
}

readMesh()

MoFEMErrorCode NonlinearPoisson::readMesh ( )

private

Definition at line 114 of file nonlinear_poisson_2d.cpp.

                                          {
  MoFEMFunctionBegin;
 
  CHKERR mField.getInterface(simpleInterface);
  CHKERR simpleInterface->getOptions();
  CHKERR simpleInterface->loadFile();
 
  MoFEMFunctionReturn(0);
}

runProgram()

MoFEMErrorCode NonlinearPoisson::runProgram

(

)

Definition at line 100 of file nonlinear_poisson_2d.cpp.

                                            {
  MoFEMFunctionBegin;
 
  readMesh();
  setupProblem();
  setIntegrationRules();
  boundaryCondition();
  assembleSystem();
  solveSystem();
  outputResults();
 
  MoFEMFunctionReturn(0);
}

setIntegrationRules()

MoFEMErrorCode NonlinearPoisson::setIntegrationRules ( )

private

Definition at line 141 of file nonlinear_poisson_2d.cpp.

                                                     {
  MoFEMFunctionBegin;
 
  auto domain_rule_lhs = [](int, int, int p) -> int { return 2 * (p + 1); };
  auto domain_rule_rhs = [](int, int, int p) -> int { return 2 * (p + 1); };
  domainTangentMatrixPipeline->getRuleHook = domain_rule_lhs;
  domainResidualVectorPipeline->getRuleHook = domain_rule_rhs;
 
  auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p + 1; };
  auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p + 1; };
  boundaryTangentMatrixPipeline->getRuleHook = boundary_rule_lhs;
  boundaryResidualVectorPipeline->getRuleHook = boundary_rule_rhs;
 
  MoFEMFunctionReturn(0);
}

setupProblem()

MoFEMErrorCode NonlinearPoisson::setupProblem ( )

private

Definition at line 124 of file nonlinear_poisson_2d.cpp.

                                              {
  MoFEMFunctionBegin;
 
  CHKERR simpleInterface->addDomainField(domainField, H1,
                                         AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
  CHKERR simpleInterface->addBoundaryField(domainField, H1,
                                           AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
 
  int order = 3;
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
  CHKERR simpleInterface->setFieldOrder(domainField, order);
 
  CHKERR simpleInterface->setUp();
 
  MoFEMFunctionReturn(0);
}

solveSystem()

MoFEMErrorCode NonlinearPoisson::solveSystem ( )

private

Definition at line 305 of file nonlinear_poisson_2d.cpp.

                                             {
  MoFEMFunctionBegin;
 
  // Create RHS and solution vectors
  SmartPetscObj<Vec> global_rhs, global_solution;
  CHKERR DMCreateGlobalVector_MoFEM(dM, global_rhs);
  global_solution = vectorDuplicate(global_rhs);
 
  // Create nonlinear solver (SNES)
  snesSolver = createSNES(mField.get_comm());
  CHKERR SNESSetFromOptions(snesSolver);
 
  // Fieldsplit block solver: yes/no
  if (1) {
    KSP ksp_solver;
    CHKERR SNESGetKSP(snesSolver, &ksp_solver);
    PC pc;
    CHKERR KSPGetPC(ksp_solver, &pc);
 
    PetscBool is_pcfs = PETSC_FALSE;
    PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
 
    // Set up FIELDSPLIT, only when option used -pc_type fieldsplit
    if (is_pcfs == PETSC_TRUE) {
      IS is_boundary;
      const MoFEM::Problem *problem_ptr;
      CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
      CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
          problem_ptr->getName(), ROW, domainField, 0, 1, &is_boundary,
          &boundaryEntitiesForFieldsplit);
      // CHKERR ISView(is_boundary, PETSC_VIEWER_STDOUT_SELF);
 
      CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
 
      CHKERR ISDestroy(&is_boundary);
    }
  }
 
  CHKERR SNESSetDM(snesSolver, dM);
  CHKERR SNESSetUp(snesSolver);
 
  // Solve the system
  CHKERR SNESSolve(snesSolver, global_rhs, global_solution);
  // VecView(global_rhs, PETSC_VIEWER_STDOUT_SELF);
 
  // Scatter result data on the mesh
  CHKERR DMoFEMMeshToGlobalVector(dM, global_solution, INSERT_VALUES,
                                  SCATTER_REVERSE);
 
  MoFEMFunctionReturn(0);
}

sourceTermFunction()

static double NonlinearPoisson::sourceTermFunction ( const double  x, const double  y, const double  z  )

inlinestaticprivate

Definition at line 28 of file nonlinear_poisson_2d.cpp.

                                                   {
    return 200 * sin(x * 10.) * cos(y * 10.);
    // return 1;
  }

Member Data Documentation

boundaryEntitiesForFieldsplit

Range NonlinearPoisson::boundaryEntitiesForFieldsplit

private

Definition at line 78 of file nonlinear_poisson_2d.cpp.

boundaryMarker

boost::shared_ptr<std::vector<unsigned char> > NonlinearPoisson::boundaryMarker

private

Definition at line 59 of file nonlinear_poisson_2d.cpp.

boundaryResidualVectorPipeline

boost::shared_ptr<EdgeEle> NonlinearPoisson::boundaryResidualVectorPipeline

private

Definition at line 65 of file nonlinear_poisson_2d.cpp.

boundaryTangentMatrixPipeline

boost::shared_ptr<EdgeEle> NonlinearPoisson::boundaryTangentMatrixPipeline

private

Definition at line 64 of file nonlinear_poisson_2d.cpp.

dM

SmartPetscObj<DM> NonlinearPoisson::dM

private

Definition at line 51 of file nonlinear_poisson_2d.cpp.

domainField

std::string NonlinearPoisson::domainField

private

Definition at line 55 of file nonlinear_poisson_2d.cpp.

domainResidualVectorPipeline

boost::shared_ptr<FaceEle> NonlinearPoisson::domainResidualVectorPipeline

private

Definition at line 63 of file nonlinear_poisson_2d.cpp.

domainTangentMatrixPipeline

boost::shared_ptr<FaceEle> NonlinearPoisson::domainTangentMatrixPipeline

private

Definition at line 62 of file nonlinear_poisson_2d.cpp.

fieldGradPtr

boost::shared_ptr<MatrixDouble> NonlinearPoisson::fieldGradPtr

private

Definition at line 70 of file nonlinear_poisson_2d.cpp.

fieldValuePtr

boost::shared_ptr<VectorDouble> NonlinearPoisson::fieldValuePtr

private

Definition at line 69 of file nonlinear_poisson_2d.cpp.

mField

MoFEM::Interface& NonlinearPoisson::mField

private

Definition at line 42 of file nonlinear_poisson_2d.cpp.

mpiComm

MPI_Comm NonlinearPoisson::mpiComm

private

Definition at line 46 of file nonlinear_poisson_2d.cpp.

mpiRank

const int NonlinearPoisson::mpiRank

private

Definition at line 48 of file nonlinear_poisson_2d.cpp.

order

int NonlinearPoisson::order

private

Definition at line 56 of file nonlinear_poisson_2d.cpp.

postProc

boost::shared_ptr<PostProcEle> NonlinearPoisson::postProc

private

Definition at line 75 of file nonlinear_poisson_2d.cpp.

previousUpdate

boost::shared_ptr<DataAtGaussPoints> NonlinearPoisson::previousUpdate

private

Definition at line 68 of file nonlinear_poisson_2d.cpp.

simpleInterface

Simple* NonlinearPoisson::simpleInterface

private

Definition at line 43 of file nonlinear_poisson_2d.cpp.

snesSolver

SmartPetscObj<SNES> NonlinearPoisson::snesSolver

private

Definition at line 52 of file nonlinear_poisson_2d.cpp.

---

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

• users_modules/tutorials/scl-4/nonlinear_poisson_2d.cpp