19 {
20
21
23
24
25 moab::Core mb_instance;
26 moab::Interface &moab = mb_instance;
27
28 try {
29
32
33 PetscInt order_x = 2;
34 PetscInt order_X = 2;
35 PetscBool flg = PETSC_TRUE;
36
37 PetscBool test_jacobian = PETSC_FALSE;
39 PETSC_NULL);
41 &flg);
43 &flg);
44
46
48
52 3);
54 3);
56
62
63 Range triangle_springs;
65 if (it->getName().compare(0, 9, "SPRING_BC") == 0) {
67 triangle_springs, true);
68 }
69 }
70
71
73 "MESH_NODE_POSITIONS");
75 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS", triangle_springs);
79
80
81 DM dm;
83
84 PetscRandom rctx;
85 PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
86
87 auto set_coord = [&](
VectorAdaptor &&field_data,
double *x,
double *y,
88 double *z) {
90 double value;
92 PetscRandomGetValue(rctx, &value);
93 field_data[0] = (*x) + (value - 0.5) *
scale;
94 PetscRandomGetValue(rctx, &value);
95 field_data[1] = (*y) + (value - 0.5) *
scale;
96 PetscRandomGetValue(rctx, &value);
97 field_data[2] = (*z) + (value - 0.5) *
scale;
99 };
100
102 "SPATIAL_POSITION");
104 set_coord, "MESH_NODE_POSITIONS");
105
106 PetscRandomDestroy(&rctx);
107
108 boost::shared_ptr<NeumannForcesSurface> surfacePressure(
110
111 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_rhs_ptr(
112 surfacePressure, &(surfacePressure->getLoopFe()));
113 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_lhs_ptr(
114 surfacePressure, &(surfacePressure->getLoopFeLhs()));
115 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_mat_rhs_ptr(
116 surfacePressure, &(surfacePressure->getLoopFeMatRhs()));
117 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_mat_lhs_ptr(
118 surfacePressure, &(surfacePressure->getLoopFeMatLhs()));
119
124
125 boost::shared_ptr<NeumannForcesSurface> surfaceForce(
127
128 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
129 fe_rhs_surface_force_ptr(surfaceForce, &(surfaceForce->getLoopFe()));
130 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
131 fe_lhs_surface_force_ptr(surfaceForce, &(surfaceForce->getLoopFeLhs()));
132 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
133 fe_mat_rhs_surface_force_ptr(surfaceForce,
134 &(surfaceForce->getLoopFeMatRhs()));
135 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
136 fe_mat_lhs_surface_force_ptr(surfaceForce,
137 &(surfaceForce->getLoopFeMatLhs()));
138
140 false, false);
142 false, false);
144 false, false);
146 false, false);
147
149 CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes,
false);
150
151 nodes.pop_front();
152 nodes.pop_back();
153
154 boost::shared_ptr<NeumannForcesSurface::DataAtIntegrationPts> dataAtPts =
155 boost::make_shared<NeumannForcesSurface::DataAtIntegrationPts>();
156
157 dataAtPts->forcesOnlyOnEntitiesRow = nodes;
158
160 if (
bit->getName().compare(0, 8,
"PRESSURE") == 0) {
161 CHKERR surfacePressure->addPressure(
"SPATIAL_POSITION", PETSC_NULL,
162 bit->getMeshsetId(),
true,
true);
163 CHKERR surfacePressure->addPressureAle(
164 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", dataAtPts,
166 true, true);
167 }
168 }
169
170 const string block_set_force_name("FORCE");
171
174 CHKERR surfaceForce->addForce(
"SPATIAL_POSITION", PETSC_NULL,
175 (
bit->getMeshsetId()),
true,
false);
176 CHKERR surfaceForce->addForceAle(
177 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", dataAtPts,
179 true, false, false);
180 }
181
183 if (it->getName().compare(0, block_set_force_name.length(),
184 block_set_force_name) == 0) {
185 CHKERR surfaceForce->addForce(
"SPATIAL_POSITION", PETSC_NULL,
186 (it->getMeshsetId()), true, true);
187 CHKERR surfaceForce->addForceAle(
188 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", dataAtPts,
190 true, true, false);
191 }
192 }
193
194
195
196 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
198 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
201 m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
202 "MESH_NODE_POSITIONS");
203
204 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ale_ptr_dx(
206 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ale_ptr_dX(
208
209 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ale_ptr(
211
212 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings>
213 data_at_spring_gp =
214 boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(
215 m_field);
216
217 Range spring_ale_nodes;
218 CHKERR moab.get_connectivity(triangle_springs, spring_ale_nodes,
true);
219
220 data_at_spring_gp->forcesOnlyOnEntitiesRow = spring_ale_nodes;
221
223 m_field, fe_spring_lhs_ale_ptr_dx, fe_spring_lhs_ale_ptr_dX,
224 fe_spring_rhs_ale_ptr, data_at_spring_gp, "SPATIAL_POSITION",
226
228 PETSC_NULL);
230 PETSC_NULL);
231
233 PETSC_NULL, PETSC_NULL);
235 PETSC_NULL, PETSC_NULL);
237 PETSC_NULL, PETSC_NULL);
238
240 nullptr, nullptr);
242 nullptr, nullptr);
243
245 nullptr, nullptr);
247 nullptr, nullptr);
248
249
251 fe_lhs_surface_force_ptr, nullptr, nullptr);
253 fe_rhs_surface_force_ptr, nullptr, nullptr);
254
256 fe_mat_lhs_surface_force_ptr, nullptr,
257 nullptr);
259 fe_mat_rhs_surface_force_ptr, nullptr,
260 nullptr);
261
263 CHKERR DMCreateGlobalVector(dm, &x);
264 CHKERR VecDuplicate(x, &f);
265 CHKERR VecSetOption(f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
267
269 CHKERR DMCreateMatrix(dm, &A);
270 CHKERR MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &fdA);
271
272 if (test_jacobian == PETSC_TRUE) {
273 char testing_options[] =
274 "-snes_test_jacobian -snes_test_jacobian_display "
275 "-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 -snes_max_it 1 ";
276
277 CHKERR PetscOptionsInsertString(NULL, testing_options);
278 } else {
279 char testing_options[] = "-snes_no_convergence_test -snes_atol 0 "
280 "-snes_rtol 0 "
281 "-snes_max_it 1 ";
282
283 CHKERR PetscOptionsInsertString(NULL, testing_options);
284 }
285
286 SNES snes;
287 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
292 CHKERR SNESSetFromOptions(snes);
293
294 CHKERR SNESSolve(snes, NULL, x);
295
296 if (test_jacobian == PETSC_FALSE) {
297 double nrm_A0;
298 CHKERR MatNorm(A, NORM_INFINITY, &nrm_A0);
299
300 char testing_options_fd[] = "-snes_fd";
301 CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
302
305 CHKERR SNESSetFromOptions(snes);
306
307 CHKERR SNESSolve(snes, NULL, x);
308 CHKERR MatAXPY(A, -1, fdA, SUBSET_NONZERO_PATTERN);
309
310 double nrm_A;
311 CHKERR MatNorm(A, NORM_INFINITY, &nrm_A);
312 PetscPrintf(PETSC_COMM_WORLD, "Matrix norms %3.4e %3.4e\n", nrm_A,
313 nrm_A / nrm_A0);
314 nrm_A /= nrm_A0;
315
316 const double tol = 1e-7;
319 "Difference between hand-calculated tangent matrix and finite "
320 "difference matrix is too big");
321 }
322 }
323
328 CHKERR SNESDestroy(&snes);
329
330
332 }
334
335
337
338 return 0;
339}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#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.
VectorShallowArrayAdaptor< double > VectorAdaptor
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.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
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.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
virtual moab::Interface & get_moab()=0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Simple interface for fast problem set-up.
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.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
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.
std::vector< std::string > & getOtherFiniteElements()
Get the Other Finite Elements.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getDomainFEName() const
Get the Domain FE Name.
Interface for nonlinear (SNES) solver.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Finite element and operators to apply force/pressures applied to surfaces.