v0.14.0
Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
MinimalSurfaceEqn Struct Reference
Collaboration diagram for MinimalSurfaceEqn:
[legend]

Public Member Functions

 MinimalSurfaceEqn (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProgram ()
 

Private Member Functions

MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setIntegrationRules ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 

Static Private Member Functions

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

Private Attributes

MoFEM::InterfacemField
 
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
 
Range boundaryEnts
 

Detailed Description

Definition at line 176 of file minimal_surface_equation.cpp.

Constructor & Destructor Documentation

◆ MinimalSurfaceEqn()

MinimalSurfaceEqn::MinimalSurfaceEqn ( MoFEM::Interface m_field)

Definition at line 207 of file minimal_surface_equation.cpp.

208  : mField(m_field) {}

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode MinimalSurfaceEqn::assembleSystem ( )
private

Definition at line 307 of file minimal_surface_equation.cpp.

307  {
309  auto add_domain_base_ops = [&](auto &pipeline) {
310  auto det_ptr = boost::make_shared<VectorDouble>();
311  auto jac_ptr = boost::make_shared<MatrixDouble>();
312  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
313  pipeline.push_back(new OpCalculateHOJac<2>(jac_ptr));
314  pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
315  pipeline.push_back(new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
316  pipeline.push_back(new OpSetHOWeightsOnFace());
317  };
318 
319  auto add_domain_lhs_ops = [&](auto &pipeline) {
320  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
321  auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
322  pipeline.push_back(
323  new OpCalculateScalarFieldGradient<2>("U", grad_u_at_gauss_pts));
324  pipeline.push_back(
325  new OpDomainTangentMatrix("U", "U", grad_u_at_gauss_pts));
326  pipeline.push_back(new OpUnSetBc("U"));
327  };
328 
329  auto add_domain_rhs_ops = [&](auto &pipeline) {
330  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
331  auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
332  pipeline.push_back(
333  new OpCalculateScalarFieldGradient<2>("U", grad_u_at_gauss_pts));
334  pipeline.push_back(new OpDomainResidualVector("U", grad_u_at_gauss_pts));
335  pipeline.push_back(new OpUnSetBc("U"));
336  };
337 
338  auto add_boundary_base_ops = [&](auto &pipeline) {};
339 
340  auto add_lhs_base_ops = [&](auto &pipeline) {
341  pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
342  pipeline.push_back(new OpBoundaryMass(
343  "U", "U", [](const double, const double, const double) { return 1; }));
344  pipeline.push_back(new OpUnSetBc("U"));
345  };
346  auto add_rhs_base_ops = [&](auto &pipeline) {
347  pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
348  auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
349  pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
350  pipeline.push_back(new OpBoundaryTimeScalarField(
351  "U", u_at_gauss_pts,
352  [](const double, const double, const double) { return 1; }));
353  pipeline.push_back(new OpBoundarySource("U", boundaryFunction));
354  pipeline.push_back(new OpUnSetBc("U"));
355  };
356 
357  auto pipeline_mng = mField.getInterface<PipelineManager>();
358  add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
359  add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
360  add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
361  add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
362 
363  add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
364  add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
365  add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
366  add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
367 
369 }

◆ boundaryCondition()

MoFEMErrorCode MinimalSurfaceEqn::boundaryCondition ( )
private

Definition at line 266 of file minimal_surface_equation.cpp.

266  {
268 
269  auto *simple = mField.getInterface<Simple>();
270 
271  auto get_ents_on_mesh_skin = [&]() {
272  Range body_ents;
273  CHKERR mField.get_moab().get_entities_by_dimension(0, 2, body_ents);
274  Skinner skin(&mField.get_moab());
275  Range skin_ents;
276  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
277  // filter not owned entities, those are not on boundary
278  Range boundary_ents;
279  ParallelComm *pcomm =
280  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
281  pcomm->filter_pstatus(skin_ents, PSTATUS_SHARED | PSTATUS_MULTISHARED,
282  PSTATUS_NOT, -1, &boundary_ents);
283 
284  Range skin_verts;
285  mField.get_moab().get_connectivity(boundary_ents, skin_verts, true);
286  boundary_ents.merge(skin_verts);
287  boundaryEnts = boundary_ents;
288  return boundary_ents;
289  };
290 
291  auto mark_boundary_dofs = [&](Range &&skin_edges) {
292  auto problem_manager = mField.getInterface<ProblemsManager>();
293  auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
294  problem_manager->markDofs(simple->getProblemName(), ROW,
295  ProblemsManager::OR, skin_edges, *marker_ptr);
296  return marker_ptr;
297  };
298 
299  // Get global local vector of marked DOFs. Is global, since is set for all
300  // DOFs on processor. Is local since only DOFs on processor are in the
301  // vector. To access DOFs use local indices.
302  boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
303 
305 }

◆ boundaryFunction()

static double MinimalSurfaceEqn::boundaryFunction ( const double  x,
const double  y,
const double  z 
)
inlinestaticprivate

Definition at line 194 of file minimal_surface_equation.cpp.

195  {
196  return sin(2 * M_PI * (x + y));
197  }

◆ outputResults()

MoFEMErrorCode MinimalSurfaceEqn::outputResults ( )
private

Definition at line 424 of file minimal_surface_equation.cpp.

424  {
426 
427  auto post_proc = boost::make_shared<PostProcEle>(mField);
428 
429  auto u_ptr = boost::make_shared<VectorDouble>();
430  post_proc->getOpPtrVector().push_back(
431  new OpCalculateScalarFieldValues("U", u_ptr));
432 
434 
435  post_proc->getOpPtrVector().push_back(
436 
437  new OpPPMap(post_proc->getPostProcMesh(),
438  post_proc->getMapGaussPts(),
439 
440  {{"U", u_ptr}},
441 
442  {}, {}, {}
443 
444  )
445 
446  );
447 
448  auto *simple = mField.getInterface<Simple>();
449  auto dm = simple->getDM();
450  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), post_proc);
451 
452  CHKERR post_proc->writeFile("out_result.h5m");
453 
455 }

◆ readMesh()

MoFEMErrorCode MinimalSurfaceEqn::readMesh ( )
private

Definition at line 224 of file minimal_surface_equation.cpp.

224  {
226 
227  auto *simple = mField.getInterface<Simple>();
229  CHKERR simple->getOptions();
230  CHKERR simple->loadFile();
231 
233 }

◆ runProgram()

MoFEMErrorCode MinimalSurfaceEqn::runProgram ( )

Definition at line 210 of file minimal_surface_equation.cpp.

210  {
212 
213  readMesh();
214  setupProblem();
217  assembleSystem();
218  solveSystem();
219  outputResults();
220 
222 }

◆ setIntegrationRules()

MoFEMErrorCode MinimalSurfaceEqn::setIntegrationRules ( )
private

Definition at line 250 of file minimal_surface_equation.cpp.

250  {
252 
253  auto integration_rule = [](int o_row, int o_col, int approx_order) {
254  return 2 * approx_order;
255  };
256 
257  auto *pipeline_mng = mField.getInterface<PipelineManager>();
258  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
259  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
260  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
261  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
262 
264 }

◆ setupProblem()

MoFEMErrorCode MinimalSurfaceEqn::setupProblem ( )
private

Definition at line 235 of file minimal_surface_equation.cpp.

235  {
237 
238  auto *simple = mField.getInterface<Simple>();
239  CHKERR simple->addDomainField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
240  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
241 
242  int order = 3;
243  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
244  CHKERR simple->setFieldOrder("U", order);
245  CHKERR simple->setUp();
246 
248 }

◆ solveSystem()

MoFEMErrorCode MinimalSurfaceEqn::solveSystem ( )
private

Definition at line 371 of file minimal_surface_equation.cpp.

371  {
373 
374  auto *simple = mField.getInterface<Simple>();
375 
376  auto set_fieldsplit_preconditioner = [&](auto snes) {
378  KSP ksp;
379  CHKERR SNESGetKSP(snes, &ksp);
380  PC pc;
381  CHKERR KSPGetPC(ksp, &pc);
382  PetscBool is_pcfs = PETSC_FALSE;
383  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
384 
385  if (is_pcfs == PETSC_TRUE) {
386  auto bc_mng = mField.getInterface<BcManager>();
387  auto name_prb = simple->getProblemName();
388  SmartPetscObj<IS> is_all_bc;
389  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
390  name_prb, ROW, "U", 0, 1, is_all_bc, &boundaryEnts);
391  int is_all_bc_size;
392  CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
393  MOFEM_LOG("EXAMPLE", Sev::inform)
394  << "Field split block size " << is_all_bc_size;
395  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
396  is_all_bc); // boundary block
397  }
399  };
400 
401  // Create global RHS and solution vectors
402  auto dm = simple->getDM();
403  SmartPetscObj<Vec> global_rhs, global_solution;
404  CHKERR DMCreateGlobalVector_MoFEM(dm, global_rhs);
405  global_solution = vectorDuplicate(global_rhs);
406 
407  // Create nonlinear solver (SNES)
408  auto pipeline_mng = mField.getInterface<PipelineManager>();
409  auto solver = pipeline_mng->createSNES();
410  CHKERR SNESSetFromOptions(solver);
411  CHKERR set_fieldsplit_preconditioner(solver);
412  CHKERR SNESSetUp(solver);
413 
414  // Solve the system
415  CHKERR SNESSolve(solver, global_rhs, global_solution);
416 
417  // Scatter result data on the mesh
418  CHKERR DMoFEMMeshToGlobalVector(dm, global_solution, INSERT_VALUES,
419  SCATTER_REVERSE);
420 
422 }

Member Data Documentation

◆ boundaryEnts

Range MinimalSurfaceEqn::boundaryEnts
private

Definition at line 204 of file minimal_surface_equation.cpp.

◆ boundaryMarker

boost::shared_ptr<std::vector<unsigned char> > MinimalSurfaceEqn::boundaryMarker
private

Definition at line 203 of file minimal_surface_equation.cpp.

◆ mField

MoFEM::Interface& MinimalSurfaceEqn::mField
private

Definition at line 200 of file minimal_surface_equation.cpp.


The documentation for this struct was generated from the following file:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MinimalSurfaceEqn::mField
MoFEM::Interface & mField
Definition: minimal_surface_equation.cpp:200
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::PipelineManager::createSNES
SmartPetscObj< SNES > createSNES(SmartPetscObj< DM > dm=nullptr)
Create SNES (nonlinear) solver.
Definition: PipelineManager.cpp:109
OpDomainTangentMatrix
Integrate the domain tangent matrix (LHS)
Definition: minimal_surface_equation.cpp:40
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
OpBoundaryTimeScalarField
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
Definition: minimal_surface_equation.cpp:16
MinimalSurfaceEqn::setupProblem
MoFEMErrorCode setupProblem()
Definition: minimal_surface_equation.cpp:235
OpDomainResidualVector
Integrate the domain residual vector (RHS)
Definition: minimal_surface_equation.cpp:118
MinimalSurfaceEqn::boundaryCondition
MoFEMErrorCode boundaryCondition()
Definition: minimal_surface_equation.cpp:266
MoFEM::OpSetHOInvJacToScalarBases< 2 >
Definition: HODataOperators.hpp:78
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MinimalSurfaceEqn::readMesh
MoFEMErrorCode readMesh()
Definition: minimal_surface_equation.cpp:224
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateHOJac< 2 >
Definition: HODataOperators.hpp:273
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
ROW
@ ROW
Definition: definitions.h:123
approx_order
static constexpr int approx_order
Definition: prism_polynomial_approximation.cpp:14
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1171
MinimalSurfaceEqn::outputResults
MoFEMErrorCode outputResults()
Definition: minimal_surface_equation.cpp:424
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MinimalSurfaceEqn::boundaryFunction
static double boundaryFunction(const double x, const double y, const double z)
Definition: minimal_surface_equation.cpp:194
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
MinimalSurfaceEqn::solveSystem
MoFEMErrorCode solveSystem()
Definition: minimal_surface_equation.cpp:371
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:539
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::OpSetHOWeightsOnFace
Modify integration weights on face to take in account higher-order geometry.
Definition: HODataOperators.hpp:122
MinimalSurfaceEqn::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: minimal_surface_equation.cpp:203
MoFEM::ProblemsManager::markDofs
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
Definition: ProblemsManager.cpp:3196
MinimalSurfaceEqn::boundaryEnts
Range boundaryEnts
Definition: minimal_surface_equation.cpp:204
MinimalSurfaceEqn::assembleSystem
MoFEMErrorCode assembleSystem()
Definition: minimal_surface_equation.cpp:307
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
MoFEM::OpUnSetBc
Definition: FormsIntegrators.hpp:49
MinimalSurfaceEqn::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
Definition: minimal_surface_equation.cpp:250
Range
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:217
OpBoundaryMass
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
Definition: minimal_surface_equation.cpp:14
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3254
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::SmartPetscObj< IS >
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:590
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
OpBoundarySource
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: minimal_surface_equation.cpp:18
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698