v0.13.2
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 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 */
29
31
32 /// \brief definition of volume element
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 */
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
77
78 // Essential boundary conditions
80
81 int oRder; ///< approximation order
82
83 // Petsc data
84 DM dM;
85 Mat A;
86 Vec D, F;
87
89 : fieldName("MAGNETIC_POTENTIAL"), feName("MAGNETIC"),
90 feNaturalBCName("MAGENTIC_NATURAL_BC"), mU(1), ePsilon(0.1) {}
92 };
93
95
96 /**
97 * \brief get natural boundary conditions
98 * \ingroup maxwell_element
99 * @return error code
100 */
101 MoFEMErrorCode getNaturalBc() {
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 }
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);
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 }
151 }
152
153 /**
154 * \brief build problem data structures
155 * \ingroup maxwell_element
156 * @return error code
157 */
158 MoFEMErrorCode createFields() {
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
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
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
203
204 // Set geometry approximation ordered
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
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
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 */
234 MoFEMErrorCode createElements() {
236 // //Elements
245 "MESH_NODE_POSITIONS");
254 blockData.feNaturalBCName, "MESH_NODE_POSITIONS");
259 // build finite elemnts
261 // build adjacencies
262 CHKERR mField.build_adjacencies(BitRefLevel().set(0));
264 }
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() {
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(
298
300 }
301
302 /**
303 * \brief destroy Distributed mesh manager
304 * \ingroup maxwell_element
305 * @return [description]
306 */
307 MoFEMErrorCode destroyProblem() {
309 CHKERR DMDestroy(&blockData.dM);
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) {
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 addHOOpsVol("MESH_NODE_POSITIONS", vol_fe, false, true, false, true);
327 vol_fe.getOpPtrVector().push_back(new OpCurlCurl(blockData));
328 vol_fe.getOpPtrVector().push_back(new OpStab(blockData));
329 TriFE tri_fe(mField);
330 tri_fe.meshPositionsFieldName = "none";
331
332 tri_fe.getOpPtrVector().push_back(
333 new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
334 tri_fe.getOpPtrVector().push_back(
335 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
336 tri_fe.getOpPtrVector().push_back(
337 new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
338 tri_fe.getOpPtrVector().push_back(new OpNaturalBC(blockData));
339
340 // create matrices and vectors
341 CHKERR DMCreateGlobalVector(blockData.dM, &blockData.D);
342 // CHKERR DMCreateMatrix(blockData.dM, &blockData.A);
343 CHKERR DMCreateGlobalVector(blockData.dM, &blockData.F);
344 CHKERR VecDuplicate(blockData.F, &blockData.D);
345
346 const MoFEM::Problem *problem_ptr;
347 CHKERR DMMoFEMGetProblemPtr(blockData.dM, &problem_ptr);
348 CHKERR mField.getInterface<MatrixManager>()
349 ->createMPIAIJ<PetscGlobalIdx_mi_tag>(problem_ptr->getName(),
350 &blockData.A);
351
352 CHKERR MatZeroEntries(blockData.A);
353 CHKERR VecZeroEntries(blockData.F);
354 CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_FORWARD);
355 CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_FORWARD);
356
357 CHKERR DMoFEMLoopFiniteElements(blockData.dM, blockData.feName.c_str(),
358 &vol_fe);
359 CHKERR DMoFEMLoopFiniteElements(blockData.dM,
360 blockData.feNaturalBCName.c_str(), &tri_fe);
361
362 CHKERR MatAssemblyBegin(blockData.A, MAT_FINAL_ASSEMBLY);
363 CHKERR MatAssemblyEnd(blockData.A, MAT_FINAL_ASSEMBLY);
364 CHKERR VecGhostUpdateBegin(blockData.F, ADD_VALUES, SCATTER_REVERSE);
365 CHKERR VecGhostUpdateEnd(blockData.F, ADD_VALUES, SCATTER_REVERSE);
366 CHKERR VecAssemblyBegin(blockData.F);
367 CHKERR VecAssemblyEnd(blockData.F);
368
369 KSP solver;
370 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
371 CHKERR KSPSetOperators(solver, blockData.A, blockData.A);
372 CHKERR KSPSetFromOptions(solver);
373 CHKERR KSPSetUp(solver);
374 CHKERR KSPSolve(solver, blockData.F, blockData.D);
375 CHKERR KSPDestroy(&solver);
376
377 CHKERR VecGhostUpdateBegin(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
378 CHKERR VecGhostUpdateEnd(blockData.D, INSERT_VALUES, SCATTER_FORWARD);
379 CHKERR DMoFEMMeshToLocalVector(blockData.dM, blockData.D, INSERT_VALUES,
380 SCATTER_REVERSE);
381
382 if (regression_test) {
383 // This test is for order = 1 only
384 double nrm2;
385 CHKERR VecNorm(blockData.D, NORM_2, &nrm2);
386 const double expected_value = 4.6772e+01;
387 const double tol = 1e-4;
388 if ((nrm2 - expected_value) / expected_value > tol) {
389 SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
390 "Regression test failed %6.4e != %6.4e", nrm2, expected_value);
391 }
392 }
393
394 CHKERR VecDestroy(&blockData.D);
395 CHKERR VecDestroy(&blockData.F);
396 CHKERR MatDestroy(&blockData.A);
397
399 }
400
401 /**
402 * \brief post-process results, i.e. save solution on the mesh
403 * \ingroup maxwell_element
404 * @return [description]
405 */
406 MoFEMErrorCode postProcessResults() {
407
409 PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
410 mField);
411
412 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, false, true, false,
413 true);
414
415 auto pos_ptr = boost::make_shared<MatrixDouble>();
416 auto field_val_ptr = boost::make_shared<MatrixDouble>();
417
418 post_proc.getOpPtrVector().push_back(
419 new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
420 post_proc.getOpPtrVector().push_back(
421 new OpCalculateHVecVectorField<3>(blockData.fieldName, field_val_ptr));
422
423 using OpPPMap = OpPostProcMapInMoab<3, 3>;
424
425 post_proc.getOpPtrVector().push_back(
426
427 new OpPPMap(
428
429 post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
430
431 {},
432
433 {{"MESH_NODE_POSITIONS", pos_ptr},
434 {blockData.fieldName, field_val_ptr}},
435
436 {},
437
438 {}
439
440 )
441
442 );
443
444 post_proc.getOpPtrVector().push_back(new OpPostProcessCurl(
445 blockData, post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
447 &post_proc);
448 CHKERR post_proc.writeFile("out_values.h5m");
450 }
451
452 /** \brief calculate and assemble CurlCurl matrix
453 * \ingroup maxwell_element
454
455 \f[
456 \mathbf{A} = \int_\Omega \mu^{-1} \left( \nabla \times \mathbf{u} \cdot
457 \nabla \times \mathbf{v} \right) \textrm{d}\Omega \f] where \f[ \mathbf{u} =
458 \nabla \times \mathbf{B} \f] where \f$\mathbf{B}\f$ is magnetic flux and
459 \f$\mu\f$ is magnetic permeability.
460
461 For more details pleas look to \cite ivanyshyn2013computation
462
463 */
466
469 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
470 data.fieldName, UserDataOperator::OPROWCOL),
471 blockData(data) {
472 sYmm = true;
473 }
474
475 MatrixDouble entityLocalMatrix;
476
477 /**
478 * \brief integrate matrix
479 * @param row_side local number of entity on element for row of the matrix
480 * @param col_side local number of entity on element for col of the matrix
481 * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
482 * @param col_type type of col entity (EDGE/TRIANGLE/TETRAHEDRON)
483 * @param row_data structure of data, like base functions and associated
484 * methods to access those data on rows
485 * @param col_data structure of data, like base functions and associated
486 * methods to access those data on rows
487 * @return error code
488 */
489 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
490 EntityType col_type,
491 EntitiesFieldData::EntData &row_data,
492 EntitiesFieldData::EntData &col_data) {
494
495 if (row_type == MBVERTEX)
497 if (col_type == MBVERTEX)
499
500 const int nb_row_dofs = row_data.getN().size2() / 3;
501 if (nb_row_dofs == 0)
503 const int nb_col_dofs = col_data.getN().size2() / 3;
504 if (nb_col_dofs == 0)
506 entityLocalMatrix.resize(nb_row_dofs, nb_col_dofs, false);
507 entityLocalMatrix.clear();
508
509 if (nb_row_dofs != static_cast<int>(row_data.getFieldData().size())) {
510 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
511 "Number of base functions and DOFs on entity is different on "
512 "rows %d!=%d",
513 nb_row_dofs, row_data.getFieldData().size());
514 }
515 if (nb_col_dofs != static_cast<int>(col_data.getFieldData().size())) {
516 SETERRQ2(
517 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
518 "Number of base functions and DOFs on entity is different on cols",
519 nb_col_dofs, col_data.getFieldData().size());
520 }
521
522 MatrixDouble row_curl_mat, col_curl_mat;
526
527 const double c0 = 1. / blockData.mU;
528 const int nb_gauss_pts = row_data.getN().size1();
529 auto t_row_curl_base = row_data.getFTensor2DiffN<3, 3>();
530
531 for (int gg = 0; gg != nb_gauss_pts; gg++) {
532
533 // get integration weight scaled by volume
534 double w = getGaussPts()(3, gg) * getVolume();
535
537 for (int aa = 0; aa != nb_row_dofs; aa++) {
538 t_row_curl(i) = levi_civita(j, i, k) * t_row_curl_base(j, k);
539
540 FTensor::Tensor0<double *> t_local_mat(&entityLocalMatrix(aa, 0), 1);
542
543 auto t_col_curl_base = col_data.getFTensor2DiffN<3, 3>(gg, 0);
544 for (int bb = 0; bb != nb_col_dofs; bb++) {
545 t_col_curl(i) = levi_civita(j, i, k) * t_col_curl_base(j, k);
546 t_local_mat += c0 * w * t_row_curl(i) * t_col_curl(i);
547 ++t_local_mat;
548 ++t_col_curl_base;
549 }
550
551 ++t_row_curl_base;
552 }
553 }
554
555 CHKERR MatSetValues(blockData.A, nb_row_dofs, &row_data.getIndices()[0],
556 nb_col_dofs, &col_data.getIndices()[0],
557 &entityLocalMatrix(0, 0), ADD_VALUES);
558
559 if (row_side != col_side || row_type != col_type) {
560 entityLocalMatrix = trans(entityLocalMatrix);
561 CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
562 nb_row_dofs, &row_data.getIndices()[0],
563 &entityLocalMatrix(0, 0), ADD_VALUES);
564 }
565
567 }
568 };
569
570 /** \brief calculate and assemble stabilization matrix
571 * \ingroup maxwell_element
572
573 \f[
574 \mathbf{A} = \int_\Omega \epsilon \mathbf{u} \cdot \mathbf{v}
575 \textrm{d}\Omega \f] where \f$\epsilon\f$ is regularization parameter.
576
577 */
578 struct OpStab
580
583 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
584 data.fieldName, UserDataOperator::OPROWCOL),
585 blockData(data) {
586 sYmm = true;
587 }
588
589 MatrixDouble entityLocalMatrix;
590
591 /**
592 * \brief integrate matrix
593 * @param row_side local number of entity on element for row of the matrix
594 * @param col_side local number of entity on element for col of the matrix
595 * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
596 * @param col_type type of col entity (EDGE/TRIANGLE/TETRAHEDRON)
597 * @param row_data structure of data, like base functions and associated
598 * methods to access those data on rows
599 * @param col_data structure of data, like base functions and associated
600 * methods to access those data on rows
601 * @return error code
602 */
603 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
604 EntityType col_type,
605 EntitiesFieldData::EntData &row_data,
606 EntitiesFieldData::EntData &col_data) {
608
609 if (row_type == MBVERTEX)
611 if (col_type == MBVERTEX)
613
614 const int nb_row_dofs = row_data.getN().size2() / 3;
615 if (nb_row_dofs == 0)
617 const int nb_col_dofs = col_data.getN().size2() / 3;
618 if (nb_col_dofs == 0)
620 entityLocalMatrix.resize(nb_row_dofs, nb_col_dofs, false);
621 entityLocalMatrix.clear();
622
623 if (nb_row_dofs != static_cast<int>(row_data.getFieldData().size())) {
624 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
625 "Number of base functions and DOFs on entity is different on "
626 "rows %d!=%d",
627 nb_row_dofs, row_data.getFieldData().size());
628 }
629 if (nb_col_dofs != static_cast<int>(col_data.getFieldData().size())) {
630 SETERRQ2(
631 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
632 "Number of base functions and DOFs on entity is different on cols",
633 nb_col_dofs, col_data.getFieldData().size());
634 }
635
636 MatrixDouble row_curl_mat, col_curl_mat;
638
639 const double c0 = 1. / blockData.mU;
640 const double c1 = blockData.ePsilon * c0;
641 const int nb_gauss_pts = row_data.getN().size1();
642
643 for (int gg = 0; gg != nb_gauss_pts; gg++) {
644
645 // get integration weight scaled by volume
646 double w = getGaussPts()(3, gg) * getVolume();
647
649 &row_data.getVectorN<3>(gg)(0, HVEC0),
650 &row_data.getVectorN<3>(gg)(0, HVEC1),
651 &row_data.getVectorN<3>(gg)(0, HVEC2), 3);
652
653 for (int aa = 0; aa != nb_row_dofs; aa++) {
654 FTensor::Tensor0<double *> t_local_mat(&entityLocalMatrix(aa, 0), 1);
656 &col_data.getVectorN<3>(gg)(0, HVEC0),
657 &col_data.getVectorN<3>(gg)(0, HVEC1),
658 &col_data.getVectorN<3>(gg)(0, HVEC2), 3);
659 for (int bb = 0; bb != nb_col_dofs; bb++) {
660 t_local_mat += c1 * w * t_row_base(i) * t_col_base(i);
661 ++t_col_base;
662 ++t_local_mat;
663 }
664 ++t_row_base;
665 }
666 }
667
668 // cerr << entityLocalMatrix << endl;
669 // cerr << endl;
670
671 CHKERR MatSetValues(blockData.A, nb_row_dofs, &row_data.getIndices()[0],
672 nb_col_dofs, &col_data.getIndices()[0],
673 &entityLocalMatrix(0, 0), ADD_VALUES);
674
675 if (row_side != col_side || row_type != col_type) {
676 entityLocalMatrix = trans(entityLocalMatrix);
677 CHKERR MatSetValues(blockData.A, nb_col_dofs, &col_data.getIndices()[0],
678 nb_row_dofs, &row_data.getIndices()[0],
679 &entityLocalMatrix(0, 0), ADD_VALUES);
680 }
681
683 }
684 };
685
686 /** \brief calculate essential boundary conditions
687 * \ingroup maxwell_element
688
689 \f[
690 \mathbf{F} = \int_{\partial\Omega} \ \mathbf{u} \cdot \mathbf{j}_i
691 \textrm{d}{\partial\Omega} \f] where \f$\mathbf{j}_i\f$ is current density
692 function.
693
694 Here simple current on coil is hard coded. In future more general
695 implementation is needed.
696
697 */
700
702
705 data.fieldName, UserDataOperator::OPROW),
706 blockData(data) {}
707
708 VectorDouble naturalBC;
709
710 /**
711 * \brief integrate matrix
712 * \ingroup maxwell_element
713 * @param row_side local number of entity on element for row of the matrix
714 * @param row_type type of row entity (EDGE/TRIANGLE/TETRAHEDRON)
715 * @param row_data structure of data, like base functions and associated
716 * methods to access those data on rows
717 * @return error code
718 */
719 MoFEMErrorCode doWork(int row_side, EntityType row_type,
720 EntitiesFieldData::EntData &row_data) {
722
723 if (row_type == MBVERTEX)
725
726 const int nb_row_dofs = row_data.getN().size2() / 3;
727 if (nb_row_dofs == 0)
729 naturalBC.resize(nb_row_dofs, false);
730 naturalBC.clear();
731
733
734 const int nb_gauss_pts = row_data.getN().size1();
735 auto t_row_base = row_data.getFTensor1N<3>();
736
737 for (int gg = 0; gg != nb_gauss_pts; gg++) {
738
739 // get integration weight scaled by volume
740 double area;
741 area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
742 double w = getGaussPts()(2, gg) * area;
743
744 // Current is on surface where natural bc are applied. It is set that
745 // current is in XY plane, circular, around the coil.
746 const double x = getCoordsAtGaussPts()(gg, 0);
747 const double y = getCoordsAtGaussPts()(gg, 1);
748 const double r = sqrt(x * x + y * y);
750 t_j(0) = -y / r;
751 t_j(1) = +x / r;
752 t_j(2) = 0;
753
754 // double a = t_j(i) * t_tangent1(i);
755 // double b = t_j(i) * t_tangent2(i);
756 // t_j(i) = a * t_tangent1(i) + b * t_tangent2(i);
757
758 // ++t_tangent1;
759 // ++t_tangent2;
760
761 FTensor::Tensor0<double *> t_f(&naturalBC[0]);
762 for (int aa = 0; aa != nb_row_dofs; aa++) {
763 t_f += w * t_row_base(i) * t_j(i);
764 ++t_row_base;
765 ++t_f;
766 }
767 }
768
769 CHKERR VecSetValues(blockData.F, row_data.getIndices().size(),
770 &row_data.getIndices()[0], &naturalBC[0], ADD_VALUES);
771
773 }
774 };
775
776 /** \brief calculate and assemble CurlCurl matrix
777 * \ingroup maxwell_element
778 */
781
783 moab::Interface &postProcMesh;
784 std::vector<EntityHandle> &mapGaussPts;
785
786 OpPostProcessCurl(BlockData &data, moab::Interface &post_proc_mesh,
787 std::vector<EntityHandle> &map_gauss_pts)
788 : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
789 data.fieldName, UserDataOperator::OPROW),
790 blockData(data), postProcMesh(post_proc_mesh),
791 mapGaussPts(map_gauss_pts) {}
792
793 MoFEMErrorCode doWork(int row_side, EntityType row_type,
794 EntitiesFieldData::EntData &row_data) {
796
797 if (row_type == MBVERTEX)
799
800 Tag th;
801 double def_val[] = {0, 0, 0};
802 CHKERR postProcMesh.tag_get_handle("MAGNETIC_INDUCTION_FIELD", 3,
803 MB_TYPE_DOUBLE, th,
804 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
805 const int nb_row_dofs = row_data.getN().size2() / 3;
806 if (nb_row_dofs == 0)
808 const void *tags_ptr[mapGaussPts.size()];
809
810 if (nb_row_dofs != row_data.getFieldData().size())
811 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
812 "Wrong number of base functions and DOFs %d != %d",
813 nb_row_dofs, row_data.getFieldData().size());
814
818 const int nb_gauss_pts = row_data.getN().size1();
819 if (nb_gauss_pts != static_cast<int>(mapGaussPts.size())) {
820 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
821 "Inconsistency number of dofs %d!=%d", nb_gauss_pts,
822 mapGaussPts.size());
823 }
824 CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0],
825 mapGaussPts.size(), tags_ptr);
826 auto t_curl_base = row_data.getFTensor2DiffN<3, 3>();
827 for (int gg = 0; gg != nb_gauss_pts; gg++) {
828 // get pointer to tag values on entity (i.e. vertex on refined
829 // post-processing mesh)
830 double *ptr = &((double *)tags_ptr[gg])[0];
832 &ptr[2]);
833 // calculate curl value
834 auto t_dof = row_data.getFTensor0FieldData();
835 for (int aa = 0; aa != nb_row_dofs; aa++) {
836 t_curl(i) += t_dof * (levi_civita(j, i, k) * t_curl_base(j, k));
837 ++t_curl_base;
838 ++t_dof;
839 }
840 }
842 }
843 };
844};
845
846#endif //__MAGNETICELEMENT_HPP__
847
848/******************************************************************************
849 * \defgroup maxwell_element Magnetic/Maxwell element
850 * \ingroup user_modules
851 ******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#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
@ HVEC0
Definition: definitions.h:186
@ HVEC1
Definition: definitions.h:186
@ HVEC2
Definition: definitions.h:186
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#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 DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:554
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.
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.
#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.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double tol
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
data structure storing material constants, model parameters, matrices, etc.
int oRder
approximation order
double mU
magnetic constant N / A2
double ePsilon
regularization paramater
calculate and assemble CurlCurl matrix
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate matrix
calculate essential boundary conditions
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
integrate matrix
calculate and assemble CurlCurl matrix
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
BlockData & blockData
OpPostProcessCurl(BlockData &data, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
calculate and assemble stabilization matrix
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate matrix
define surface element
TriFE(MoFEM::Interface &m_field)
definition of volume element
VolumeFE(MoFEM::Interface &m_field)
Implementation of magneto-static problem (basic Implementation)
MoFEMErrorCode createFields()
build problem data structures
MoFEMErrorCode createElements()
create finite elements
MoFEMErrorCode destroyProblem()
destroy Distributed mesh manager
MoFEMErrorCode postProcessResults()
post-process results, i.e. save solution on the mesh
MoFEMErrorCode createProblem()
create problem
MoFEM::Interface & mField
MoFEMErrorCode getEssentialBc()
get essential boundary conditions (only homogenous case is considered)
MoFEMErrorCode getNaturalBc()
get natural boundary conditions
MagneticElement(MoFEM::Interface &m_field)
virtual ~MagneticElement()=default
MoFEMErrorCode solveProblem(const bool regression_test=false)
solve problem
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.
Deprecated interface functions.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Post post-proc data at points from hash maps.
keeps basic data about problem
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)