v0.14.0
Loading...
Searching...
No Matches
Functions | Variables
minimal_surface_area.cpp File Reference
#include <MoFEM.hpp>
#include <MethodForForceScaling.hpp>
#include <DirichletBC.hpp>
#include <PostProcOnRefMesh.hpp>
#include <MinimalSurfaceElement.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 41 of file minimal_surface_area.cpp.

41 {
42
43 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
44
45 try {
46
47 moab::Core mb_instance;
48 moab::Interface &moab = mb_instance;
49
50 char mesh_file_name[255];
51 PetscInt order = 2;
52 PetscInt nb_sub_steps = 10;
53 PetscBool flg_file = PETSC_FALSE;
54 PetscBool is_partitioned = PETSC_FALSE;
55 PetscBool is_atom_test = PETSC_FALSE;
56
57 {
58
59 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "",
60 "Minimal surface area config", "none");
61 CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
62 mesh_file_name, 255, &flg_file);
63 CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
64 &order, PETSC_NULL);
65 CHKERR PetscOptionsInt("-my_nb_sub_steps", "number of substeps", "", 10,
66 &nb_sub_steps, PETSC_NULL);
67 CHKERR PetscOptionsBool("-my_is_partitioned",
68 "set if mesh is partitioned (this result that "
69 "each process keeps only part of the mesh",
70 "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
71 CHKERR PetscOptionsBool(
72 "-my_is_atom_test",
73 "is used with testing, exit with error when diverged", "",
74 PETSC_FALSE, &is_atom_test, PETSC_NULL);
75
76 ierr = PetscOptionsEnd();
78
79 // Reade parameters from line command
80 if (flg_file != PETSC_TRUE) {
81 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
82 }
83
84 if (is_partitioned == PETSC_TRUE) {
85 // Read mesh to MOAB
86 const char *option;
87 option = "PARALLEL=BCAST_DELETE;"
88 "PARALLEL_RESOLVE_SHARED_ENTS;"
89 "PARTITION=PARALLEL_PARTITION;";
90 CHKERR moab.load_file(mesh_file_name, 0, option);
91 } else {
92 const char *option;
93 option = "";
94 CHKERR moab.load_file(mesh_file_name, 0, option);
95 }
96 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
97 if (pcomm == NULL)
98 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
99 }
100
101 MoFEM::Core core(moab);
102 MoFEM::Interface &m_field = core;
103 // meshset consisting all entities in mesh
104 EntityHandle root_set = moab.get_root_set();
105
106 // Seed all mesh entities to MoFEM database, those entities can be
107 // potentially used as finite elements or as entities which carry some
108 // approximation field.
109 BitRefLevel bit_level0;
110 bit_level0.set(0);
111 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
112 0, 2, bit_level0);
113
114 // Add field
115 CHKERR m_field.add_field("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
116 // Add entities to field (root_mesh, i.e. on all mesh entities fields are
117 // approx.) On those entities and lower dimension entities approximation is
118 // spaned,
119 CHKERR m_field.add_ents_to_field_by_type(root_set, MBTRI, "U");
120 // Set approximation oder
121 CHKERR m_field.set_field_order(root_set, MBTRI, "U", order);
122 CHKERR m_field.set_field_order(root_set, MBEDGE, "U", order);
123 CHKERR m_field.set_field_order(root_set, MBVERTEX, "U", 1);
124 CHKERR m_field.build_fields();
125
126 // Add finite element (this defines element declaration comes later)
127 CHKERR m_field.add_finite_element("MINIMAL_SURFACE_ELEMENT");
128 // Set on which fields that element operates
130 "MINIMAL_SURFACE_ELEMENT", "U");
132 "MINIMAL_SURFACE_ELEMENT", "U");
134 "MINIMAL_SURFACE_ELEMENT", "U");
135 // Add entities to that element
137 root_set, MBTRI, "MINIMAL_SURFACE_ELEMENT");
138
139 CHKERR m_field.add_finite_element("BC_ELEMENT");
140 CHKERR m_field.modify_finite_element_add_field_row("BC_ELEMENT", "U");
141 CHKERR m_field.modify_finite_element_add_field_col("BC_ELEMENT", "U");
142 CHKERR m_field.modify_finite_element_add_field_data("BC_ELEMENT", "U");
143 // Add entities to that element
144 Range bc_edges, bc_ents;
145 {
146 // Take a skin form the mesh
147 Skinner skin(&m_field.get_moab());
148 Range tris;
149 CHKERR moab.get_entities_by_type(root_set, MBTRI, tris, false);
150 Range bc_edges;
151 CHKERR skin.find_skin(root_set, tris, false, bc_edges);
152 Range bc_nodes;
153 CHKERR moab.get_connectivity(bc_edges, bc_nodes);
154 bc_ents.merge(bc_edges);
155 bc_ents.merge(bc_nodes);
156 // Add entities to that element
157 CHKERR m_field.add_ents_to_finite_element_by_type(bc_edges, MBEDGE,
158 "BC_ELEMENT");
159 }
160
161 // build finite elements
163 // build adjacencies between elements and degrees of freedom
164 CHKERR m_field.build_adjacencies(bit_level0);
165
166 // Create minimal surface area and elements
167 MinimalSurfaceElement minimal_surface_element(m_field);
168 // Data structure for common data
169 MinimalSurfaceElement::CommonData mse_common_data;
170
171 // This class is used to fix dofs on entities
172 DirichletFixFieldAtEntitiesBc fix_edges_ents(m_field, "U", NULL, NULL, NULL,
173 bc_ents);
174
175 CHKERR DMRegister_MoFEM("MoFEM");
176 // Solve BC problem
177 {
178 // define problems
179 DM bc_dm;
180 // Register DM problem
181 CHKERR DMCreate(PETSC_COMM_WORLD, &bc_dm);
182 CHKERR DMSetType(bc_dm, "MoFEM");
183 // Create DM instance
184 CHKERR DMMoFEMCreateMoFEM(bc_dm, &m_field, "DM_BC", bit_level0);
185 // Configure DM form line command options (DM itself, solvers,
186 // pre-conditioners, ... )
187 CHKERR DMSetFromOptions(bc_dm);
188 // Add elements to dm (only one here)
189 CHKERR DMMoFEMAddElement(bc_dm, "BC_ELEMENT");
190 // Set up problem (number dofs, partition mesh, etc.)
191 CHKERR DMSetUp(bc_dm);
192 Mat A;
193 Vec T, F;
194 CHKERR DMCreateGlobalVector(bc_dm, &T);
195 CHKERR VecDuplicate(T, &F);
196 CHKERR DMCreateMatrix(bc_dm, &A);
197
198 // Boundary problem
199 minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
201 minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
203 CHKERR DMoFEMLoopFiniteElements(bc_dm, "BC_ELEMENT",
204 &minimal_surface_element.feBcEdge);
205 CHKERR VecAssemblyBegin(F);
206 CHKERR VecAssemblyEnd(F);
207 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
208 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
209
210 {
211 KSP solver;
212 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
213 CHKERR KSPSetOperators(solver, A, A);
214 CHKERR KSPSetFromOptions(solver);
215 CHKERR KSPSetUp(solver);
216 CHKERR KSPSolve(solver, F, T);
217 CHKERR KSPDestroy(&solver);
218 }
219
220 // Scatter data on mesh
221 CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
222 CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
223 CHKERR DMoFEMMeshToGlobalVector(bc_dm, T, INSERT_VALUES, SCATTER_REVERSE);
224
225 CHKERR VecDestroy(&T);
226 CHKERR VecDestroy(&F);
227 CHKERR MatDestroy(&A);
228 CHKERR DMDestroy(&bc_dm);
229 }
230
231 // //define problems
232 // CHKERR m_field.add_problem("DM_MINIMAL_AREA");
233 // //set refinement level for problem
234 // CHKERR
235 // m_field.modify_problem_ref_level_add_bit("DM_MINIMAL_AREA",bit_level0);
236 // Create discrete manager instance
237 DM dm;
238 {
239 // Register problem
240 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
241 CHKERR DMSetType(dm, "MoFEM");
242 // Create DM instance
243 CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "DM_MINIMAL_AREA", bit_level0);
244 // Get options fron line command
245 CHKERR DMSetFromOptions(dm);
246 // Add elements to dm (only one here)
247 CHKERR DMMoFEMAddElement(dm, "MINIMAL_SURFACE_ELEMENT");
248 // Set up data structures
249 CHKERR DMSetUp(dm);
250 }
251
252 // Create matrices and vectors used for analysis
253 Vec T, F;
254 Mat A;
255 {
256 CHKERR DMCreateGlobalVector(dm, &T);
257 CHKERR VecDuplicate(T, &F);
258 CHKERR DMCreateMatrix(dm, &A);
259 }
260
261 // Adding elements to DMSnes
262 {
263
264 // Right hand side operators (run for each entity in element)
265 // Calculate inverse of jacobian
266 auto det_ptr = boost::make_shared<VectorDouble>();
267 auto jac_ptr = boost::make_shared<MatrixDouble>();
268 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
269 minimal_surface_element.feRhs.getOpPtrVector().push_back(
270 new OpCalculateHOJacForFace(jac_ptr));
271 minimal_surface_element.feRhs.getOpPtrVector().push_back(
272 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
273 // Apply inv_jac_ptr to shape functions
274 minimal_surface_element.feRhs.getOpPtrVector().push_back(
275 new OpSetInvJacH1ForFace(inv_jac_ptr));
276 // Below operators are specific to minimal area problem
277 minimal_surface_element.feRhs.getOpPtrVector().push_back(
278 new MinimalSurfaceElement::OpGetDataAtGaussPts("U", mse_common_data));
279 minimal_surface_element.feRhs.getOpPtrVector().push_back(
281 "U", mse_common_data, true));
282 minimal_surface_element.feRhs.getOpPtrVector().push_back(
283 new MinimalSurfaceElement::OpAssembleResidaul("U", mse_common_data));
284
285 // Left hand side operators (run for each entity in element)
286 // Calculate inverse of jacobian
287 minimal_surface_element.feLhs.getOpPtrVector().push_back(
288 new OpCalculateHOJacForFace(jac_ptr));
289 minimal_surface_element.feLhs.getOpPtrVector().push_back(
290 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
291 // Apply inv_jac_ptr to shape functions
292 minimal_surface_element.feLhs.getOpPtrVector().push_back(
293 new OpSetInvJacH1ForFace(inv_jac_ptr));
294 // Below operators are specific to minimal area problem
295 minimal_surface_element.feLhs.getOpPtrVector().push_back(
296 new MinimalSurfaceElement::OpGetDataAtGaussPts("U", mse_common_data));
297 minimal_surface_element.feLhs.getOpPtrVector().push_back(
299 "U", mse_common_data, false));
300 minimal_surface_element.feLhs.getOpPtrVector().push_back(
301 new MinimalSurfaceElement::OpAssembleTangent("U", mse_common_data));
302
303 // Rhs
304 CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL, &fix_edges_ents,
305 NULL);
306 CHKERR DMMoFEMSNESSetFunction(dm, "MINIMAL_SURFACE_ELEMENT",
307 &minimal_surface_element.feRhs, PETSC_NULL,
308 PETSC_NULL);
310 &fix_edges_ents);
311
312 // Lhs
313 CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, NULL, &fix_edges_ents,
314 NULL);
315 CHKERR DMMoFEMSNESSetJacobian(dm, "MINIMAL_SURFACE_ELEMENT",
316 &minimal_surface_element.feLhs, PETSC_NULL,
317 PETSC_NULL);
319 &fix_edges_ents);
320 }
321
322 SNESConvergedReason snes_reason;
323 // Solve problem
324 {
325 SNES snes;
326 SnesCtx *snes_ctx;
327 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
328 // CHKERR SNESSetDM(snes,dm);
329 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
330 CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
331 CHKERR SNESSetJacobian(snes, A, A, SnesMat, snes_ctx);
332 CHKERR SNESSetFromOptions(snes);
333
334 Vec T0;
335 CHKERR VecDuplicate(T, &T0);
336 CHKERR DMoFEMMeshToLocalVector(dm, T0, INSERT_VALUES, SCATTER_FORWARD);
337
338 double step_size = 1. / nb_sub_steps;
339 for (int ss = 1; ss <= nb_sub_steps; ss++) {
340 CHKERR VecAXPY(T, step_size, T0);
341 CHKERR DMoFEMMeshToLocalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
342 CHKERR SNESSolve(snes, PETSC_NULL, T);
343 CHKERR SNESGetConvergedReason(snes, &snes_reason);
344 int its;
345 CHKERR SNESGetIterationNumber(snes, &its);
346 CHKERR PetscPrintf(PETSC_COMM_WORLD,
347 "number of Newton iterations = %D\n\n", its);
348 if (snes_reason < 0) {
349 if (is_atom_test) {
350 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
351 "atom test diverged");
352 } else {
353 break;
354 }
355 }
356 }
357 CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
358 CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
359 CHKERR DMoFEMMeshToGlobalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
360
361 // CHKERR MatView(A,PETSC_VIEWER_DRAW_SELF);
362 // std::string wait;
363 // std::cin >> wait;
364
365 CHKERR VecDestroy(&T0);
366 CHKERR SNESDestroy(&snes);
367 }
368
369 {
370 PostProcFaceOnRefinedMesh post_proc(m_field);
371 CHKERR post_proc.generateReferenceElementMesh();
372 CHKERR post_proc.addFieldValuesPostProc("U");
373 CHKERR DMoFEMLoopFiniteElements(dm, "MINIMAL_SURFACE_ELEMENT",
374 &post_proc);
375 // Save data on mesh
376 CHKERR post_proc.postProcMesh.write_file("out.h5m", "MOAB",
377 "PARALLEL=WRITE_PART");
378 }
379
380 // Clean and destroy
381 {
382 CHKERR DMDestroy(&dm);
383 CHKERR VecDestroy(&T);
384 CHKERR VecDestroy(&F);
385 CHKERR MatDestroy(&A);
386 }
387 }
389
391 return 0;
392}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
@ F
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
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 DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
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 DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1080
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:521
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
char mesh_file_name[255]
static char help[]
const double T
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:136
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:27
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
constexpr AssemblyType A
Fix dofs on entities.
Evaluate function values and gradients at Gauss Pts.
Implementation of minimal area element.
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
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.
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Postprocess on face.

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 39 of file minimal_surface_area.cpp.