v0.16.3
Loading...
Searching...
No Matches
MagneticElement.hpp
Go to the documentation of this file.
1/** \file MagneticElement.hpp
2 * \brief Implementation of magnetic element
3 * \ingroup maxwell_element
4 * \example mofem/tutorials/max-0_magnetostatics/src/MagneticElement.hpp
5 *
6 */
7
8#ifndef __MAGNETICELEMENT_HPP__
9#define __MAGNETICELEMENT_HPP__
10
11/**
12 * \brief Implementation of magneto-static problem (basic Implementation)
13 * \ingroup maxwell_element
14 *
15 * Look for theory and details here:
16 *
17 * \cite ivanyshyn2013computation
18 * <www.hpfem.jku.at/publications/szthesis.pdf>
19 * <https://pdfs.semanticscholar.org/0574/e5763d9b64bff16f908e9621f23af3b3dc86.pdf>
20 *
21 * Election file and all other problem related file are here \ref
22 * maxwell_element.
23 *
24 * \todo Extension for mix formulation
25 * \todo Use appropriate pre-conditioner for large problems
26 *
27 */
28struct MagneticElement {
29
30 MoFEM::Interface &mField;
31
32 /// \brief definition of volume element
33 struct VolumeFE : public MoFEM::VolumeElementForcesAndSourcesCore {
36 int getRule(int order) { return 2 * order + 1; };
37 };
38
39 // /// \brief definition of volume element
40 // struct VolumeFEReducedIntegration: public
41 // MoFEM::VolumeElementForcesAndSourcesCore {
42 // VolumeFEReducedIntegration(MoFEM::Interface &m_field):
43 // MoFEM::VolumeElementForcesAndSourcesCore(m_field) {}
44 // int getRule(int order) { return 2*order+1; };
45 // };
46
47 /** \brief define surface element
48 *
49 */
50 struct TriFE : public MoFEM::FaceElementForcesAndSourcesCore {
53 int getRule(int order) { return 2 * order; };
54 };
55
56 MagneticElement(MoFEM::Interface &m_field) : mField(m_field) {}
57 virtual ~MagneticElement() = default;
58
59 /**
60 * \brief data structure storing material constants, model parameters,
61 * matrices, etc.
62 *
63 */
64 struct BlockData {
65
66 // field
67 const string fieldName;
68 const string feName;
69 const string feNaturalBCName;
70
71 // material parameters
72 double mU; ///< magnetic constant N / A2
73 double ePsilon; ///< regularization paramater
74
75 // Natural boundary conditions
76 Range naturalBc;
77
78 // Essential boundary conditions
79 Range essentialBc;
80
81 int oRder; ///< approximation order
82
83 // Petsc data
84 DM dM;
85 Mat A;
86 Vec D, F;
87
88 BlockData()
89 : fieldName("MAGNETIC_POTENTIAL"), feName("MAGNETIC"),
90 feNaturalBCName("MAGENTIC_NATURAL_BC"), mU(1), ePsilon(0.1) {}
91 ~BlockData() {}
92 };
93
94 BlockData blockData;
95
96 /**
97 * \brief get natural boundary conditions
98 * \ingroup maxwell_element
99 * @return error code
100 */
101 MoFEMErrorCode getNaturalBc() {
102 MoFEMFunctionBegin;
103 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, bit)) {
104 if (bit->getName().compare(0, 9, "NATURALBC") == 0) {
105 Range faces;
106 CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI,
107 faces, true);
108 CHKERR mField.get_moab().get_adjacencies(
109 faces, 1, true, blockData.naturalBc, moab::Interface::UNION);
110 blockData.naturalBc.merge(faces);
111 }
112 }
113 MoFEMFunctionReturn(0);
114 }
115
116 /**
117 * \brief get essential boundary conditions (only homogenous case is
118 * considered) \ingroup maxwell_element
119 * @return error code
120 */
121 MoFEMErrorCode getEssentialBc() {
122 ParallelComm *pcomm =
123 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
124 MoFEMFunctionBegin;
125 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, bit)) {
126 if (bit->getName().compare(0, 10, "ESSENTIALBC") == 0) {
127 Range faces;
128 CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI,
129 faces, true);
130 CHKERR mField.get_moab().get_adjacencies(
131 faces, 1, true, blockData.essentialBc, moab::Interface::UNION);
132 blockData.essentialBc.merge(faces);
133 }
134 }
135 if (blockData.essentialBc.empty()) {
136 Range tets;
137 CHKERR mField.get_moab().get_entities_by_type(0, MBTET, tets);
138 Skinner skin(&mField.get_moab());
139 Range skin_faces; // skin faces from 3d ents
140 CHKERR skin.find_skin(0, tets, false, skin_faces);
141 skin_faces = subtract(skin_faces, blockData.naturalBc);
142 Range proc_skin;
143 CHKERR pcomm->filter_pstatus(skin_faces,
144 PSTATUS_SHARED | PSTATUS_MULTISHARED,
145 PSTATUS_NOT, -1, &proc_skin);
146 CHKERR mField.get_moab().get_adjacencies(
147 proc_skin, 1, true, blockData.essentialBc, moab::Interface::UNION);
148 blockData.essentialBc.merge(proc_skin);
149 }
150 MoFEMFunctionReturn(0);
151 }
152
153 /**
154 * \brief build problem data structures
155 * \ingroup maxwell_element
156 * @return error code
157 */
158 MoFEMErrorCode createFields() {
159 MoFEMFunctionBegin;
160
161 // Set entities bit level. each entity has bit level depending for example
162 // on refinement level. In this case we do not refine mesh or not do
163 // topological changes, simply set refinement level to zero on all entities.
164
165 CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
166 0, 3, BitRefLevel().set(0));
167
168 // add fields
169 CHKERR mField.add_field(blockData.fieldName, HCURL, DEMKOWICZ_JACOBI_BASE,
170 1);
171 CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
172 3);
173 // meshset consisting all entities in mesh
174 EntityHandle root_set = mField.get_moab().get_root_set();
175 // add entities to field
176 CHKERR mField.add_ents_to_field_by_type(root_set, MBTET,
177 blockData.fieldName);
178
179 // // The higher-order gradients can be gauged by locally skipping the
180 // // corresponding degrees of freedom and basis functions in the
181 // higher-order
182 // // edge-based, face-based and cell-based finite element subspaces.
183 //
184 // Range tris,edges;
185 // CHKERR mField.get_moab().get_entities_by_type(root_set,MBTRI,tris,true);
186 // mField.get_moab().get_entities_by_type(root_set,MBEDGE,edges,true);
187 //
188 // // Set order in volume
189 // Range bc_ents = unite(blockData.naturalBc,blockData.essentialBc);
190 // Range vol_ents = subtract(unite(tris,edges),bc_ents);
191 // CHKERR
192 // mField.set_field_order(vol_ents,blockData.fieldName,blockData.oRder); int
193 // gauged_order = 1; CHKERR
194 // mField.set_field_order(bc_ents,blockData.fieldName,gauged_order);
195
196 // Set order on tets
197 CHKERR mField.set_field_order(root_set, MBTET, blockData.fieldName,
198 blockData.oRder);
199 CHKERR mField.set_field_order(root_set, MBTRI, blockData.fieldName,
200 blockData.oRder);
201 CHKERR mField.set_field_order(root_set, MBEDGE, blockData.fieldName,
202 blockData.oRder);
203
204 // Set geometry approximation ordered
205 CHKERR mField.add_ents_to_field_by_type(root_set, MBTET,
206 "MESH_NODE_POSITIONS");
207 CHKERR mField.set_field_order(root_set, MBTET, "MESH_NODE_POSITIONS", 2);
208 CHKERR mField.set_field_order(root_set, MBTRI, "MESH_NODE_POSITIONS", 2);
209 CHKERR mField.set_field_order(root_set, MBEDGE, "MESH_NODE_POSITIONS", 2);
210 CHKERR mField.set_field_order(root_set, MBVERTEX, "MESH_NODE_POSITIONS", 1);
211
212 // build field
213 CHKERR mField.build_fields();
214
215 // get HO geometry for 10 node tets
216 // This method takes coordinates form edges mid nodes in 10 node tet and
217 // project values on 2nd order hierarchical basis used to approx. geometry.
218 Projection10NodeCoordsOnField ent_method_material(mField,
219 "MESH_NODE_POSITIONS");
220 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
221
222 MoFEMFunctionReturn(0);
223 }
224
225 /**
226 * \brief create finite elements
227 * \ingroup maxwell_element
228 *
229 * Create volume and surface element. Surface element is used to integrate
230 * natural boundary conditions.
231 *
232 * @return error code
233 */
265
266 /**
267 * \brief create problem
268 * \ingroup maxwell_element
269 *
270 * Problem is collection of finite elements. With the information on which
271 * fields finite elements operates the matrix and left and right hand side
272 * vector could be created.
273 *
274 * Here we use Distributed mesh manager from PETSc as a inteface.
275 *
276 * @return error code
277 */
278 MoFEMErrorCode createProblem() {
279 MoFEMFunctionBegin;
280 // set up DM
281 CHKERR DMRegister_MoFEM("DMMOFEM");
282 CHKERR DMCreate(PETSC_COMM_WORLD, &blockData.dM);
283 CHKERR DMSetType(blockData.dM, "DMMOFEM");
284 CHKERR DMMoFEMCreateMoFEM(blockData.dM, &mField, "MAGNETIC_PROBLEM",
285 BitRefLevel().set(0));
286 CHKERR DMSetFromOptions(blockData.dM);
287 CHKERR DMMoFEMSetIsPartitioned(blockData.dM, PETSC_TRUE);
288 // add elements to blockData.dM
289 CHKERR DMMoFEMAddElement(blockData.dM, blockData.feName);
290 CHKERR DMMoFEMAddElement(blockData.dM, blockData.feNaturalBCName);
291 CHKERR DMSetUp(blockData.dM);
292
293 // remove essential DOFs
294 const MoFEM::Problem *problem_ptr;
295 CHKERR DMMoFEMGetProblemPtr(blockData.dM, &problem_ptr);
296 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
297 problem_ptr->getName(), blockData.fieldName, blockData.essentialBc);
298
299 MoFEMFunctionReturn(0);
300 }
301
302 /**
303 * \brief destroy Distributed mesh manager
304 * \ingroup maxwell_element
305 * @return [description]
306 */
307 MoFEMErrorCode destroyProblem() {
308 MoFEMFunctionBegin;
309 CHKERR DMDestroy(&blockData.dM);
310 MoFEMFunctionReturn(0);
311 }
312
313 /** \brief solve problem
314 * \ingroup maxwell_element
315 *
316 * Create matrices; integrate over elements; solve linear system of equations
317 *
318 */
319 MoFEMErrorCode solveProblem(const bool regression_test = false) {
320 MoFEMFunctionBegin;
321
322 VolumeFE vol_fe(mField);
323 auto material_grad_mat = boost::make_shared<MatrixDouble>();
324 auto material_det_vec = boost::make_shared<VectorDouble>();
325 auto material_inv_grad_mat = boost::make_shared<MatrixDouble>();
326 CHKERR AddHOOps<3, 3, 3>::add(vol_fe.getOpPtrVector(), {HCURL, L2},
327 "MESH_NODE_POSITIONS");
328 vol_fe.getOpPtrVector().push_back(new OpCurlCurl(blockData));
329 vol_fe.getOpPtrVector().push_back(new OpStab(blockData));
330 TriFE tri_fe(mField);
331 tri_fe.meshPositionsFieldName = "none";
332
333 tri_fe.getOpPtrVector().push_back(
334 new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
335 tri_fe.getOpPtrVector().push_back(
336 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
337 tri_fe.getOpPtrVector().push_back(
338 new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
339 tri_fe.getOpPtrVector().push_back(new OpNaturalBC(blockData));
340
341 // create matrices and vectors
342 CHKERR DMCreateGlobalVector(blockData.dM, &blockData.D);
343 // CHKERR DMCreateMatrix(blockData.dM, &blockData.A);
344 CHKERR DMCreateGlobalVector(blockData.dM, &blockData.F);
345 CHKERR VecDuplicate(blockData.F, &blockData.D);
346
347 const MoFEM::Problem *problem_ptr;
348 CHKERR DMMoFEMGetProblemPtr(blockData.dM, &problem_ptr);
349 CHKERR mField.getInterface<MatrixManager>()
350 ->createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_ptr->getName(),
351 &blockData.A);
352
353 CHKERR MatZeroEntries(blockData.A);
354 CHKERR VecZeroEntries(blockData.F);
355 CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_FORWARD);
356 CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_FORWARD);
357
358 CHKERR DMoFEMLoopFiniteElements(blockData.dM, blockData.feName.c_str(),
359 &vol_fe);
360 CHKERR DMoFEMLoopFiniteElements(blockData.dM,
361 blockData.feNaturalBCName.c_str(), &tri_fe);
362
363 CHKERR MatAssemblyBegin(blockData.A, MAT_FINAL_ASSEMBLY);
364 CHKERR MatAssemblyEnd(blockData.A, MAT_FINAL_ASSEMBLY);
365 CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_REVERSE);
366 CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_REVERSE);
367 CHKERR VecAssemblyBegin(blockData.F);
368 CHKERR VecAssemblyEnd(blockData.F);
369
370 KSP solver;
371 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
372 CHKERR KSPSetOperators(solver, blockData.A, blockData.A);
373 CHKERR KSPSetFromOptions(solver);
374 CHKERR KSPSetUp(solver);
375 CHKERR KSPSolve(solver, blockData.F, blockData.D);
376 CHKERR mField.getInterface<TestingInterface>()->checkSolver(
377 solver, "magnetostatics");
378 CHKERR KSPDestroy(&solver);
379
380 CHKERR VecGhostUpdateBegin(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
381 CHKERR VecGhostUpdateEnd(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
382 CHKERR DMoFEMMeshToLocalVector(blockData.dM, blockData.D, INSERT_VALUES,
383 SCATTER_REVERSE);
384
385 if (regression_test) {
386 // This test is for order = 1 only
387 double nrm2;
388 CHKERR VecNorm(blockData.D, NORM_2, &nrm2);
389 const double expected_value = 4.6772e+01;
390 const double tol = 1e-4;
391 if ((nrm2 - expected_value) / expected_value > tol) {
392 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
393 "Regression test failed %6.4e != %6.4e", nrm2, expected_value);
394 }
395 }
396
397 CHKERR VecDestroy(&blockData.D);
398 CHKERR VecDestroy(&blockData.F);
399 CHKERR MatDestroy(&blockData.A);
400
401 MoFEMFunctionReturn(0);
402 }
403
404 /**
405 * \brief post-process results, i.e. save solution on the mesh
406 * \ingroup maxwell_element
407 * @return [description]
408 */
409 MoFEMErrorCode postProcessResults() {
410
411 MoFEMFunctionBegin;
412 PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
413 mField);
414
415 CHKERR AddHOOps<3, 3, 3>::add(post_proc.getOpPtrVector(), {HCURL, L2},
416 "MESH_NODE_POSITIONS");
417
418 auto pos_ptr = boost::make_shared<MatrixDouble>();
419 auto field_val_ptr = boost::make_shared<MatrixDouble>();
420
421 post_proc.getOpPtrVector().push_back(
422 new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
423 post_proc.getOpPtrVector().push_back(
424 new OpCalculateHVecVectorField<3>(blockData.fieldName, field_val_ptr));
425
426 using OpPPMap = OpPostProcMapInMoab<3, 3>;
427
428 post_proc.getOpPtrVector().push_back(
429
430 new OpPPMap(
431
432 post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
433
434 {},
435
436 {{"MESH_NODE_POSITIONS", pos_ptr},
437 {blockData.fieldName, field_val_ptr}},
438
439 {},
440
441 {}
442
443 )
444
445 );
446
447 post_proc.getOpPtrVector().push_back(new OpPostProcessCurl(
448 blockData, post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
449 CHKERR DMoFEMLoopFiniteElements(blockData.dM, blockData.feName.c_str(),
450 &post_proc);
451 CHKERR post_proc.writeFile("out_values.h5m");
452 MoFEMFunctionReturn(0);
453 }
454
455 /** \brief calculate and assemble CurlCurl matrix
456 * \ingroup maxwell_element
457
458 \f[
459 \mathbf{A} = \int_\Omega \mu^{-1} \left( \nabla \times \mathbf{u} \cdot
460 \nabla \times \mathbf{v} \right) \textrm{d}\Omega \f] where \f[ \mathbf{u} =
461 \nabla \times \mathbf{B} \f] where \f$\mathbf{B}\f$ is magnetic flux and
462 \f$\mu\f$ is magnetic permeability.
463
464 For more details pleas look to \cite ivanyshyn2013computation
465
466 */
467 struct OpCurlCurl
468 : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
469
470 BlockData &blockData;
471 OpCurlCurl(BlockData &data)
472 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
473 data.fieldName, UserDataOperator::OPROWCOL),
474 blockData(data) {
475 sYmm = true;
476 }
477
478 MatrixDouble entityLocalMatrix;
479
480 /**
481 * \brief integrate matrix
482 * @param row_side local number of entity on element for row of the matrix
483 * @param col_side local number of entity on element for col of the matrix
484 * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
485 * @param col_type type of col entity (EDGE/TRIANGLE/TETRAHEDRON)
486 * @param row_data structure of data, like base functions and associated
487 * methods to access those data on rows
488 * @param col_data structure of data, like base functions and associated
489 * methods to access those data on rows
490 * @return error code
491 */
492 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
493 EntityType col_type,
494 EntitiesFieldData::EntData &row_data,
495 EntitiesFieldData::EntData &col_data) {
496 MoFEMFunctionBegin;
497
498 if (row_type == MBVERTEX)
499 MoFEMFunctionReturnHot(0);
500 if (col_type == MBVERTEX)
501 MoFEMFunctionReturnHot(0);
502
503 const int nb_row_dofs = row_data.getN().size2() / 3;
504 if (nb_row_dofs == 0)
505 MoFEMFunctionReturnHot(0);
506 const int nb_col_dofs = col_data.getN().size2() / 3;
507 if (nb_col_dofs == 0)
508 MoFEMFunctionReturnHot(0);
509 entityLocalMatrix.resize(nb_row_dofs, nb_col_dofs, false);
510 entityLocalMatrix.clear();
511
512 if (nb_row_dofs != static_cast<int>(row_data.getFieldData().size())) {
513 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
514 "Number of base functions and DOFs on entity is different on "
515 "rows %d!=%zu",
516 nb_row_dofs, row_data.getFieldData().size());
517 }
518 if (nb_col_dofs != static_cast<int>(col_data.getFieldData().size())) {
519 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
520 "Number of base functions and DOFs on entity is different on "
521 "cols %d!=%zu",
522 nb_col_dofs, col_data.getFieldData().size());
523 }
524
525 MatrixDouble row_curl_mat, col_curl_mat;
526 FTensor::Index<'i', 3> i;
527 FTensor::Index<'j', 3> j;
528 FTensor::Index<'k', 3> k;
529
530 const double c0 = 1. / blockData.mU;
531 const int nb_gauss_pts = row_data.getN().size1();
532 auto t_row_curl_base = row_data.getFTensor2DiffN<3, 3>();
533
534 for (int gg = 0; gg != nb_gauss_pts; gg++) {
535
536 // get integration weight scaled by volume
537 double w = getGaussPts()(3, gg) * getVolume();
538
539 FTensor::Tensor1<double, 3> t_row_curl;
540 for (int aa = 0; aa != nb_row_dofs; aa++) {
541 t_row_curl(i) = levi_civita(j, i, k) * t_row_curl_base(j, k);
542
543 FTensor::Tensor0<double *> t_local_mat(&entityLocalMatrix(aa, 0), 1);
544 FTensor::Tensor1<double, 3> t_col_curl;
545
546 auto t_col_curl_base = col_data.getFTensor2DiffN<3, 3>(gg, 0);
547 for (int bb = 0; bb != nb_col_dofs; bb++) {
548 t_col_curl(i) = levi_civita(j, i, k) * t_col_curl_base(j, k);
549 t_local_mat += c0 * w * t_row_curl(i) * t_col_curl(i);
550 ++t_local_mat;
551 ++t_col_curl_base;
552 }
553
554 ++t_row_curl_base;
555 }
556 }
557
558 CHKERR MatSetValues(blockData.A, nb_row_dofs, &row_data.getIndices()[0],
559 nb_col_dofs, &col_data.getIndices()[0],
560 &entityLocalMatrix(0, 0), ADD_VALUES);
561
562 if (row_side != col_side || row_type != col_type) {
563 entityLocalMatrix = trans(entityLocalMatrix);
564 CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
565 nb_row_dofs, &row_data.getIndices()[0],
566 &entityLocalMatrix(0, 0), ADD_VALUES);
567 }
568
569 MoFEMFunctionReturn(0);
570 }
571 };
572
573 /** \brief calculate and assemble stabilization matrix
574 * \ingroup maxwell_element
575
576 \f[
577 \mathbf{A} = \int_\Omega \epsilon \mathbf{u} \cdot \mathbf{v}
578 \textrm{d}\Omega \f] where \f$\epsilon\f$ is regularization parameter.
579
580 */
581 struct OpStab
582 : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
583
584 BlockData &blockData;
585 OpStab(BlockData &data)
586 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
587 data.fieldName, UserDataOperator::OPROWCOL),
588 blockData(data) {
589 sYmm = true;
590 }
591
592 MatrixDouble entityLocalMatrix;
593
594 /**
595 * \brief integrate matrix
596 * @param row_side local number of entity on element for row of the matrix
597 * @param col_side local number of entity on element for col of the matrix
598 * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
599 * @param col_type type of col entity (EDGE/TRIANGLE/TETRAHEDRON)
600 * @param row_data structure of data, like base functions and associated
601 * methods to access those data on rows
602 * @param col_data structure of data, like base functions and associated
603 * methods to access those data on rows
604 * @return error code
605 */
606 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
607 EntityType col_type,
608 EntitiesFieldData::EntData &row_data,
609 EntitiesFieldData::EntData &col_data) {
610 MoFEMFunctionBegin;
611
612 if (row_type == MBVERTEX)
613 MoFEMFunctionReturnHot(0);
614 if (col_type == MBVERTEX)
615 MoFEMFunctionReturnHot(0);
616
617 const int nb_row_dofs = row_data.getN().size2() / 3;
618 if (nb_row_dofs == 0)
619 MoFEMFunctionReturnHot(0);
620 const int nb_col_dofs = col_data.getN().size2() / 3;
621 if (nb_col_dofs == 0)
622 MoFEMFunctionReturnHot(0);
623 entityLocalMatrix.resize(nb_row_dofs, nb_col_dofs, false);
624 entityLocalMatrix.clear();
625
626 if (nb_row_dofs != static_cast<int>(row_data.getFieldData().size())) {
627 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
628 "Number of base functions and DOFs on entity is different on "
629 "rows %d!=%zu",
630 nb_row_dofs, row_data.getFieldData().size());
631 }
632 if (nb_col_dofs != static_cast<int>(col_data.getFieldData().size())) {
633 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
634 "Number of base functions and DOFs on entity is different on "
635 "cols %d!=%zu",
636 nb_col_dofs, col_data.getFieldData().size());
637 }
638
639 MatrixDouble row_curl_mat, col_curl_mat;
640 FTensor::Index<'i', 3> i;
641
642 const double c0 = 1. / blockData.mU;
643 const double c1 = blockData.ePsilon * c0;
644 const int nb_gauss_pts = row_data.getN().size1();
645
646 for (int gg = 0; gg != nb_gauss_pts; gg++) {
647
648 // get integration weight scaled by volume
649 double w = getGaussPts()(3, gg) * getVolume();
650
651 FTensor::Tensor1<const double *, 3> t_row_base(
652 &row_data.getVectorN<3>(gg)(0, HVEC0),
653 &row_data.getVectorN<3>(gg)(0, HVEC1),
654 &row_data.getVectorN<3>(gg)(0, HVEC2), 3);
655
656 for (int aa = 0; aa != nb_row_dofs; aa++) {
657 FTensor::Tensor0<double *> t_local_mat(&entityLocalMatrix(aa, 0), 1);
658 FTensor::Tensor1<const double *, 3> t_col_base(
659 &col_data.getVectorN<3>(gg)(0, HVEC0),
660 &col_data.getVectorN<3>(gg)(0, HVEC1),
661 &col_data.getVectorN<3>(gg)(0, HVEC2), 3);
662 for (int bb = 0; bb != nb_col_dofs; bb++) {
663 t_local_mat += c1 * w * t_row_base(i) * t_col_base(i);
664 ++t_col_base;
665 ++t_local_mat;
666 }
667 ++t_row_base;
668 }
669 }
670
671 // cerr << entityLocalMatrix << endl;
672 // cerr << endl;
673
674 CHKERR MatSetValues(blockData.A, nb_row_dofs, &row_data.getIndices()[0],
675 nb_col_dofs, &col_data.getIndices()[0],
676 &entityLocalMatrix(0, 0), ADD_VALUES);
677
678 if (row_side != col_side || row_type != col_type) {
679 entityLocalMatrix = trans(entityLocalMatrix);
680 CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
681 nb_row_dofs, &row_data.getIndices()[0],
682 &entityLocalMatrix(0, 0), ADD_VALUES);
683 }
684
685 MoFEMFunctionReturn(0);
686 }
687 };
688
689 /** \brief calculate essential boundary conditions
690 * \ingroup maxwell_element
691
692 \f[
693 \mathbf{F} = \int_{\partial\Omega} \ \mathbf{u} \cdot \mathbf{j}_i
694 \textrm{d}{\partial\Omega} \f] where \f$\mathbf{j}_i\f$ is current density
695 function.
696
697 Here simple current on coil is hard coded. In future more general
698 implementation is needed.
699
700 */
701 struct OpNaturalBC
702 : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
703
704 BlockData &blockData;
705
706 OpNaturalBC(BlockData &data)
707 : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
708 data.fieldName, UserDataOperator::OPROW),
709 blockData(data) {}
710
711 VectorDouble naturalBC;
712
713 /**
714 * \brief integrate matrix
715 * \ingroup maxwell_element
716 * @param row_side local number of entity on element for row of the matrix
717 * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
718 * @param row_data structure of data, like base functions and associated
719 * methods to access those data on rows
720 * @return error code
721 */
722 MoFEMErrorCode doWork(int row_side, EntityType row_type,
723 EntitiesFieldData::EntData &row_data) {
724 MoFEMFunctionBegin;
725
726 if (row_type == MBVERTEX)
727 MoFEMFunctionReturnHot(0);
728
729 const int nb_row_dofs = row_data.getN().size2() / 3;
730 if (nb_row_dofs == 0)
731 MoFEMFunctionReturnHot(0);
732 naturalBC.resize(nb_row_dofs, false);
733 naturalBC.clear();
734
735 FTensor::Index<'i', 3> i;
736
737 const int nb_gauss_pts = row_data.getN().size1();
738 auto t_row_base = row_data.getFTensor1N<3>();
739
740 for (int gg = 0; gg != nb_gauss_pts; gg++) {
741
742 // get integration weight scaled by volume
743 double area;
744 area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
745 double w = getGaussPts()(2, gg) * area;
746
747 // Current is on surface where natural bc are applied. It is set that
748 // current is in XY plane, circular, around the coil.
749 const double x = getCoordsAtGaussPts()(gg, 0);
750 const double y = getCoordsAtGaussPts()(gg, 1);
751 const double r = sqrt(x * x + y * y);
752 FTensor::Tensor1<double, 3> t_j;
753 t_j(0) = -y / r;
754 t_j(1) = +x / r;
755 t_j(2) = 0;
756
757 // double a = t_j(i) * t_tangent1(i);
758 // double b = t_j(i) * t_tangent2(i);
759 // t_j(i) = a * t_tangent1(i) + b * t_tangent2(i);
760
761 // ++t_tangent1;
762 // ++t_tangent2;
763
764 FTensor::Tensor0<double *> t_f(&naturalBC[0]);
765 for (int aa = 0; aa != nb_row_dofs; aa++) {
766 t_f += w * t_row_base(i) * t_j(i);
767 ++t_row_base;
768 ++t_f;
769 }
770 }
771
772 CHKERR VecSetValues(blockData.F, row_data.getIndices().size(),
773 &row_data.getIndices()[0], &naturalBC[0], ADD_VALUES);
774
775 MoFEMFunctionReturn(0);
776 }
777 };
778
779 /** \brief calculate and assemble CurlCurl matrix
780 * \ingroup maxwell_element
781 */
782 struct OpPostProcessCurl
783 : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
784
785 BlockData &blockData;
786 moab::Interface &postProcMesh;
787 std::vector<EntityHandle> &mapGaussPts;
788
789 OpPostProcessCurl(BlockData &data, moab::Interface &post_proc_mesh,
790 std::vector<EntityHandle> &map_gauss_pts)
791 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
792 data.fieldName, UserDataOperator::OPROW),
793 blockData(data), postProcMesh(post_proc_mesh),
794 mapGaussPts(map_gauss_pts) {}
795
796 MoFEMErrorCode doWork(int row_side, EntityType row_type,
797 EntitiesFieldData::EntData &row_data) {
798 MoFEMFunctionBegin;
799
800 if (row_type == MBVERTEX)
801 MoFEMFunctionReturnHot(0);
802
803 Tag th;
804 double def_val[] = {0, 0, 0};
805 CHKERR postProcMesh.tag_get_handle("MAGNETIC_INDUCTION_FIELD", 3,
806 MB_TYPE_DOUBLE, th,
807 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
808 const int nb_row_dofs = row_data.getN().size2() / 3;
809 if (nb_row_dofs == 0)
810 MoFEMFunctionReturnHot(0);
811 const void *tags_ptr[mapGaussPts.size()];
812
813 if (nb_row_dofs != row_data.getFieldData().size())
814 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
815 "Wrong number of base functions and DOFs %d != %zu",
816 nb_row_dofs, row_data.getFieldData().size());
817
818 FTensor::Index<'i', 3> i;
819 FTensor::Index<'j', 3> j;
820 FTensor::Index<'k', 3> k;
821 const int nb_gauss_pts = row_data.getN().size1();
822 if (nb_gauss_pts != static_cast<int>(mapGaussPts.size())) {
823 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
824 "Inconsistency number of dofs %d!=%zu", nb_gauss_pts,
825 mapGaussPts.size());
826 }
827 CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0],
828 mapGaussPts.size(), tags_ptr);
829 auto t_curl_base = row_data.getFTensor2DiffN<3, 3>();
830 for (int gg = 0; gg != nb_gauss_pts; gg++) {
831 // get pointer to tag values on entity (i.e. vertex on refined
832 // post-processing mesh)
833 double *ptr = &((double *)tags_ptr[gg])[0];
834 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_curl(ptr, &ptr[1],
835 &ptr[2]);
836 // calculate curl value
837 auto t_dof = row_data.getFTensor0FieldData();
838 for (int aa = 0; aa != nb_row_dofs; aa++) {
839 t_curl(i) += t_dof * (levi_civita(j, i, k) * t_curl_base(j, k));
840 ++t_curl_base;
841 ++t_dof;
842 }
843 }
844 MoFEMFunctionReturn(0);
845 }
846 };
847};
848
849#endif //__MAGNETICELEMENT_HPP__
850
851/******************************************************************************
852 * \defgroup maxwell_element Magnetic/Maxwell element
853 * \ingroup user_modules
854 ******************************************************************************/
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition HookeElement.hpp:75
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
H1
@ H1
continuous field
Definition definitions.h:85
HCURL
@ HCURL
field with continuous tangents
Definition definitions.h:86
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
BLOCKSET
@ BLOCKSET
Definition definitions.h:161
HVEC0
@ HVEC0
Definition definitions.h:199
HVEC1
@ HVEC1
Definition definitions.h:199
HVEC2
@ HVEC2
Definition definitions.h:199
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
order
constexpr int order
Definition dg_projection.cpp:38
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:605
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::set_field_order
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.
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition MeshsetsManager.hpp:94
bit
auto bit
set bit
Definition hanging_node_approx.cpp:75
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
tol
double tol
Definition mesh_smoothing.cpp:23
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition photon_diffusion.cpp:40
MagneticElement::BlockData
data structure storing material constants, model parameters, matrices, etc.
Definition MagneticElement.hpp:64
MagneticElement::BlockData::oRder
int oRder
approximation order
Definition MagneticElement.hpp:81
MagneticElement::BlockData::mU
double mU
magnetic constant N / A2
Definition MagneticElement.hpp:72
MagneticElement::BlockData::ePsilon
double ePsilon
regularization paramater
Definition MagneticElement.hpp:73
MagneticElement::OpCurlCurl
calculate and assemble CurlCurl matrix
Definition MagneticElement.hpp:468
MagneticElement::OpCurlCurl::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate matrix
Definition MagneticElement.hpp:492
MagneticElement::OpNaturalBC
calculate essential boundary conditions
Definition MagneticElement.hpp:702
MagneticElement::OpNaturalBC::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
integrate matrix
Definition MagneticElement.hpp:722
MagneticElement::OpPostProcessCurl
calculate and assemble CurlCurl matrix
Definition MagneticElement.hpp:783
MagneticElement::OpPostProcessCurl::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition MagneticElement.hpp:787
MagneticElement::OpPostProcessCurl::OpPostProcessCurl
OpPostProcessCurl(BlockData &data, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts)
Definition MagneticElement.hpp:789
MagneticElement::OpPostProcessCurl::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition MagneticElement.hpp:796
MagneticElement::OpStab
calculate and assemble stabilization matrix
Definition MagneticElement.hpp:582
MagneticElement::OpStab::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate matrix
Definition MagneticElement.hpp:606
MagneticElement::TriFE
define surface element
Definition MagneticElement.hpp:50
MagneticElement::TriFE::TriFE
TriFE(MoFEM::Interface &m_field)
Definition MagneticElement.hpp:51
MagneticElement::VolumeFE
definition of volume element
Definition MagneticElement.hpp:33
MagneticElement::VolumeFE::VolumeFE
VolumeFE(MoFEM::Interface &m_field)
Definition MagneticElement.hpp:34
MagneticElement
Implementation of magneto-static problem (basic Implementation)
Definition MagneticElement.hpp:28
MagneticElement::createFields
MoFEMErrorCode createFields()
build problem data structures
Definition MagneticElement.hpp:158
MagneticElement::createElements
MoFEMErrorCode createElements()
create finite elements
Definition MagneticElement.hpp:234
MagneticElement::destroyProblem
MoFEMErrorCode destroyProblem()
destroy Distributed mesh manager
Definition MagneticElement.hpp:307
MagneticElement::postProcessResults
MoFEMErrorCode postProcessResults()
post-process results, i.e. save solution on the mesh
Definition MagneticElement.hpp:409
MagneticElement::createProblem
MoFEMErrorCode createProblem()
create problem
Definition MagneticElement.hpp:278
MagneticElement::mField
MoFEM::Interface & mField
Definition MagneticElement.hpp:30
MagneticElement::getEssentialBc
MoFEMErrorCode getEssentialBc()
get essential boundary conditions (only homogenous case is considered)
Definition MagneticElement.hpp:121
MagneticElement::getNaturalBc
MoFEMErrorCode getNaturalBc()
get natural boundary conditions
Definition MagneticElement.hpp:101
MagneticElement::MagneticElement
MagneticElement(MoFEM::Interface &m_field)
Definition MagneticElement.hpp:56
MagneticElement::~MagneticElement
virtual ~MagneticElement()=default
MagneticElement::solveProblem
MoFEMErrorCode solveProblem(const bool regression_test=false)
solve problem
Definition MagneticElement.hpp:319
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::add_field
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.
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition ForcesAndSourcesCore.hpp:121
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition PostProcBrokenMeshInMoabBase.hpp:745
MoFEM::Problem
keeps basic data about problem
Definition ProblemsMultiIndices.hpp:54
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:103
Generated by Doxygen 1.9.8 and hosted at University of Glasgow