v0.14.0
Loading...
Searching...
No Matches
nonlinear_poisson_2d.cpp
Go to the documentation of this file.
1#include <stdlib.h>
4
5using namespace MoFEM;
6using namespace NonlinearPoissonOps;
7
8static char help[] = "...\n\n";
9
11public:
13
14 // Declaration of the main function to run analysis
16
17private:
18 // Declaration of other main functions called in runProgram()
26
27 // Function to calculate the Source term
28 static double sourceTermFunction(const double x, const double y,
29 const double z) {
30 return 200 * sin(x * 10.) * cos(y * 10.);
31 // return 1;
32 }
33
34 // Function to calculate the Boundary term
35 static double boundaryFunction(const double x, const double y,
36 const double z) {
37 return sin(x * 10.) * cos(y * 10.);
38 // return 0;
39 }
40
41 // Main interfaces
44
45 // mpi parallel communicator
46 MPI_Comm mpiComm;
47 // Number of processors
48 const int mpiRank;
49
50 // Discrete Manager and nonliner SNES solver using SmartPetscObj
53
54 // Field name and approximation order
55 std::string domainField;
56 int order;
57
58 // Object to mark boundary entities for the assembling of domain elements
59 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
60
61 // MoFEM working Pipelines for LHS and RHS of domain and boundary
62 boost::shared_ptr<FaceEle> domainTangentMatrixPipeline;
63 boost::shared_ptr<FaceEle> domainResidualVectorPipeline;
64 boost::shared_ptr<EdgeEle> boundaryTangentMatrixPipeline;
65 boost::shared_ptr<EdgeEle> boundaryResidualVectorPipeline;
66
67 // Objects needed for solution updates in Newton's method
68 boost::shared_ptr<DataAtGaussPoints> previousUpdate;
69 boost::shared_ptr<VectorDouble> fieldValuePtr;
70 boost::shared_ptr<MatrixDouble> fieldGradPtr;
71
73
74 // Object needed for postprocessing
75 boost::shared_ptr<PostProcEle> postProc;
76
77 // Boundary entities marked for fieldsplit (block) solver - optional
79};
80
82 : domainField("U"), mField(m_field), mpiComm(mField.get_comm()),
83 mpiRank(mField.get_comm_rank()) {
84 domainTangentMatrixPipeline = boost::shared_ptr<FaceEle>(new FaceEle(mField));
86 boost::shared_ptr<FaceEle>(new FaceEle(mField));
88 boost::shared_ptr<EdgeEle>(new EdgeEle(mField));
90 boost::shared_ptr<EdgeEle>(new EdgeEle(mField));
91
93 boost::shared_ptr<DataAtGaussPoints>(new DataAtGaussPoints());
94 fieldValuePtr = boost::shared_ptr<VectorDouble>(previousUpdate,
95 &previousUpdate->fieldValue);
96 fieldGradPtr = boost::shared_ptr<MatrixDouble>(previousUpdate,
97 &previousUpdate->fieldGrad);
98}
99
102
103 readMesh();
104 setupProblem();
108 solveSystem();
110
112}
113
116
120
122}
123
126
131
132 int order = 3;
133 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
135
137
139}
140
143
144 auto domain_rule_lhs = [](int, int, int p) -> int { return 2 * (p + 1); };
145 auto domain_rule_rhs = [](int, int, int p) -> int { return 2 * (p + 1); };
146 domainTangentMatrixPipeline->getRuleHook = domain_rule_lhs;
147 domainResidualVectorPipeline->getRuleHook = domain_rule_rhs;
148
149 auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p + 1; };
150 auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p + 1; };
151 boundaryTangentMatrixPipeline->getRuleHook = boundary_rule_lhs;
152 boundaryResidualVectorPipeline->getRuleHook = boundary_rule_rhs;
153
155}
156
159
160 // Get boundary edges marked in block named "BOUNDARY_CONDITION"
161 auto get_ents_on_mesh_skin = [&]() {
162 Range boundary_entities;
164 std::string entity_name = it->getName();
165 if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
166 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
167 boundary_entities, true);
168 }
169 }
170 // Add vertices to boundary entities
171 Range boundary_vertices;
172 CHKERR mField.get_moab().get_connectivity(boundary_entities,
173 boundary_vertices, true);
174 boundary_entities.merge(boundary_vertices);
175
176 // Store entities for fieldsplit (block) solver
177 boundaryEntitiesForFieldsplit = boundary_entities;
178
179 return boundary_entities;
180 };
181
182 auto mark_boundary_dofs = [&](Range &&skin_edges) {
183 auto problem_manager = mField.getInterface<ProblemsManager>();
184 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
185 problem_manager->markDofs(simpleInterface->getProblemName(), ROW,
186 skin_edges, *marker_ptr);
187 return marker_ptr;
188 };
189
190 // Get global local vector of marked DOFs. Is global, since is set for all
191 // DOFs on processor. Is local since only DOFs on processor are in the
192 // vector. To access DOFs use local indices.
193 boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
194
196}
197
200
201 auto det_ptr = boost::make_shared<VectorDouble>();
202 auto jac_ptr = boost::make_shared<MatrixDouble>();
203 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
204
205 { // Push operators to the Pipeline that is responsible for calculating the
206 // domain tangent matrix (LHS)
207
208 // Add default operators to calculate inverse of Jacobian (needed for
209 // implementation of 2D problem but not 3D ones)
210
211 domainTangentMatrixPipeline->getOpPtrVector().push_back(
212 new OpCalculateHOJac<2>(jac_ptr));
213 domainTangentMatrixPipeline->getOpPtrVector().push_back(
214 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
215 domainTangentMatrixPipeline->getOpPtrVector().push_back(
216 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
217 domainTangentMatrixPipeline->getOpPtrVector().push_back(
219
220 // Add default operator to calculate field values at integration points
221 domainTangentMatrixPipeline->getOpPtrVector().push_back(
223 // Add default operator to calculate field gradient at integration points
224 domainTangentMatrixPipeline->getOpPtrVector().push_back(
226
227 // Push operators for domain tangent matrix (LHS)
228 domainTangentMatrixPipeline->getOpPtrVector().push_back(
231 }
232
233 { // Push operators to the Pipeline that is responsible for calculating the
234 // domain residual vector (RHS)
235
236 // Add default operators to calculate inverse of Jacobian (needed for
237 // implementation of 2D problem but not 3D ones)
238 domainResidualVectorPipeline->getOpPtrVector().push_back(
239 new OpCalculateHOJac<2>(jac_ptr));
240 domainResidualVectorPipeline->getOpPtrVector().push_back(
241 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
242 domainResidualVectorPipeline->getOpPtrVector().push_back(
243 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
244 domainResidualVectorPipeline->getOpPtrVector().push_back(
246
247 // Add default operator to calculate field values at integration points
248 domainResidualVectorPipeline->getOpPtrVector().push_back(
250 // Add default operator to calculate field gradient at integration points
251 domainResidualVectorPipeline->getOpPtrVector().push_back(
253
254 domainResidualVectorPipeline->getOpPtrVector().push_back(
257 }
258
259 { // Push operators to the Pipeline that is responsible for calculating the
260 // boundary tangent matrix (LHS)
261
262 boundaryTangentMatrixPipeline->getOpPtrVector().push_back(
264 }
265
266 { // Push operators to the Pipeline that is responsible for calculating
267 // boundary residual vector (RHS)
268
269 // Add default operator to calculate field values at integration points
270 boundaryResidualVectorPipeline->getOpPtrVector().push_back(
272
273 boundaryResidualVectorPipeline->getOpPtrVector().push_back(
276 }
277
278 // get Discrete Manager (SmartPetscObj)
280
281 { // Set operators for nonlinear equations solver (SNES) from MoFEM Pipelines
282
283 // Set operators for calculation of LHS and RHS of domain elements
284 boost::shared_ptr<FaceEle> null_face;
287 null_face);
290 null_face);
291
292 // Set operators for calculation of LHS and RHS of boundary elements
293 boost::shared_ptr<EdgeEle> null_edge;
296 null_edge);
299 null_edge);
300 }
301
303}
304
307
308 // Create RHS and solution vectors
309 SmartPetscObj<Vec> global_rhs, global_solution;
311 global_solution = vectorDuplicate(global_rhs);
312
313 // Create nonlinear solver (SNES)
315 CHKERR SNESSetFromOptions(snesSolver);
316
317 // Fieldsplit block solver: yes/no
318 if (1) {
319 KSP ksp_solver;
320 CHKERR SNESGetKSP(snesSolver, &ksp_solver);
321 PC pc;
322 CHKERR KSPGetPC(ksp_solver, &pc);
323
324 PetscBool is_pcfs = PETSC_FALSE;
325 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
326
327 // Set up FIELDSPLIT, only when option used -pc_type fieldsplit
328 if (is_pcfs == PETSC_TRUE) {
329 IS is_boundary;
330 const MoFEM::Problem *problem_ptr;
331 CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
332 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
333 problem_ptr->getName(), ROW, domainField, 0, 1, &is_boundary,
335 // CHKERR ISView(is_boundary, PETSC_VIEWER_STDOUT_SELF);
336
337 CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
338
339 CHKERR ISDestroy(&is_boundary);
340 }
341 }
342
343 CHKERR SNESSetDM(snesSolver, dM);
344 CHKERR SNESSetUp(snesSolver);
345
346 // Solve the system
347 CHKERR SNESSolve(snesSolver, global_rhs, global_solution);
348 // VecView(global_rhs, PETSC_VIEWER_STDOUT_SELF);
349
350 // Scatter result data on the mesh
351 CHKERR DMoFEMMeshToGlobalVector(dM, global_solution, INSERT_VALUES,
352 SCATTER_REVERSE);
353
355}
356
359
360 postProc = boost::make_shared<PostProcEle>(mField);
361
362 auto u_ptr = boost::make_shared<VectorDouble>();
363 postProc->getOpPtrVector().push_back(
365
367
368 postProc->getOpPtrVector().push_back(
369
370 new OpPPMap(postProc->getPostProcMesh(), postProc->getMapGaussPts(),
371
372 {{"U", u_ptr}},
373
374 {}, {}, {}
375
376 )
377
378 );
379
381 dM, simpleInterface->getDomainFEName(),
382 boost::dynamic_pointer_cast<FEMethod>(postProc));
383
384 CHKERR postProc->writeFile("out_result.h5m");
385
387}
388
389int main(int argc, char *argv[]) {
390
391 // Initialisation of MoFEM/PETSc and MOAB data structures
392 const char param_file[] = "param_file.petsc";
394
395 // Error handling
396 try {
397 // Register MoFEM discrete manager in PETSc
398 DMType dm_name = "DMMOFEM";
399 CHKERR DMRegister_MoFEM(dm_name);
400
401 // Create MOAB instance
402 moab::Core mb_instance; // mesh database
403 moab::Interface &moab = mb_instance; // mesh database interface
404
405 // Create MoFEM instance
406 MoFEM::Core core(moab); // finite element database
407 MoFEM::Interface &m_field = core; // finite element interface
408
409 // Run the main analysis
410 NonlinearPoisson poisson_problem(m_field);
411 CHKERR poisson_problem.runProgram();
412 }
414
415 // Finish work: cleaning memory, getting statistics, etc.
417
418 return 0;
419}
static Index< 'p', 3 > p
std::string param_file
int main()
Definition: adol-c_atom.cpp:46
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
@ 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
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMoFEM.cpp:704
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:412
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMoFEM.cpp:745
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
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
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
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createSNES(MPI_Comm comm)
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
static char help[]
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
keeps basic data about problem
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
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:348
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEMErrorCode addBoundaryField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on boundary.
Definition: Simple.cpp:282
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:362
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:341
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode outputResults()
MoFEMErrorCode boundaryCondition()
static double boundaryFunction(const double x, const double y, const double z)
MoFEMErrorCode solveSystem()
SmartPetscObj< SNES > snesSolver
static double sourceTermFunction(const double x, const double y, const double z)
MoFEMErrorCode readMesh()
MoFEM::Interface & mField
SmartPetscObj< DM > dM
MoFEMErrorCode assembleSystem()
NonlinearPoisson(MoFEM::Interface &m_field)
boost::shared_ptr< EdgeEle > boundaryTangentMatrixPipeline
boost::shared_ptr< VectorDouble > fieldValuePtr
MoFEMErrorCode setIntegrationRules()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
boost::shared_ptr< MatrixDouble > fieldGradPtr
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode runProgram()
boost::shared_ptr< EdgeEle > boundaryResidualVectorPipeline
boost::shared_ptr< FaceEle > domainResidualVectorPipeline
MoFEMErrorCode setupProblem()
boost::shared_ptr< DataAtGaussPoints > previousUpdate
boost::shared_ptr< FaceEle > domainTangentMatrixPipeline