v0.13.2
Loading...
Searching...
No Matches
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}
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
static double boundaryFunction(const double x, const double y, const double z)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Integrate the domain residual vector (RHS)
Integrate the domain tangent matrix (LHS)

◆ 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}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
@ ROW
Definition: definitions.h:123
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define CHKERR
Inline error check.
Definition: definitions.h:535
virtual moab::Interface & get_moab()=0
Problem manager is used to build and partition problems.
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.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27

◆ 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}
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:572
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Post post-proc data at points from hash maps.

◆ readMesh()

MoFEMErrorCode MinimalSurfaceEqn::readMesh ( )
private

Definition at line 224 of file minimal_surface_equation.cpp.

224 {
226
227 auto *simple = mField.getInterface<Simple>();
228 CHKERR 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();
218 solveSystem();
220
222}
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode setIntegrationRules()

◆ 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}
auto integration_rule
static constexpr int approx_order

◆ 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}
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)

◆ 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}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1153
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:521
SmartPetscObj< SNES > createSNES(SmartPetscObj< DM > dm=nullptr)
Create SNES (nonlinear) solver.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
intrusive_ptr for managing petsc objects

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: