v0.14.0
Loading...
Searching...
No Matches
DirichletBC.cpp
Go to the documentation of this file.
1
2
3#include <MoFEM.hpp>
4
5using namespace MoFEM;
7#include <DirichletBC.hpp>
8
9using namespace boost::numeric;
10
12 const Problem *problem_ptr, const FieldBitNumber bit_number, Range &ents,
13 boost::function<
14 MoFEMErrorCode(const boost::shared_ptr<MoFEM::NumeredDofEntity> &dof)>
15 for_each_dof) {
17
18 auto &dofs_by_uid = problem_ptr->getNumeredRowDofsPtr()->get<Unique_mi_tag>();
19
20 for (auto eit = ents.pair_begin(); eit != ents.pair_end(); ++eit) {
21
22 auto lo_dit = dofs_by_uid.lower_bound(
23 DofEntity::getLoFieldEntityUId(bit_number, eit->first));
24 auto hi_dit = dofs_by_uid.upper_bound(
25 DofEntity::getHiFieldEntityUId(bit_number, eit->second));
26
27 for (; lo_dit != hi_dit; ++lo_dit) {
28 auto &dof = *lo_dit;
29 if (dof->getHasLocalIndex())
30 CHKERR for_each_dof(dof);
31 }
32 }
33
35};
36
38 const std::string &field_name,
39 Mat Aij, Vec X, Vec F,
40 string blockset_name)
41 : mField(m_field), fieldName(field_name), blocksetName(blockset_name),
42 dIag(1) {
43 snes_B = Aij;
44 snes_x = X;
45 snes_f = F;
46 ts_B = Aij;
47 ts_u = X;
48 ts_F = F;
49};
50
52 const std::string &field_name,
53 string blockset_name)
54 : mField(m_field), fieldName(field_name), blocksetName(blockset_name),
55 dIag(1) {
56 snes_B = PETSC_NULL;
57 snes_x = PETSC_NULL;
58 snes_f = PETSC_NULL;
59 ts_B = PETSC_NULL;
60 ts_u = PETSC_NULL;
61 ts_F = PETSC_NULL;
62};
63
65 std::vector<DataFromBc> &bc_data) {
66
68
69 // Loop over meshsets with Dirichlet boundary condition on displacements
72 bc_data.push_back(DataFromBc());
74 CHKERR bc_data.back().getBcData(mydata, &(*it));
75 CHKERR bc_data.back().getEntitiesFromBc(mField, &(*it));
76 }
77 // Loop over blocksets with DISPLACEMENT (default) boundary condition
79 if (it->getName().compare(0, blocksetName.length(), blocksetName) == 0) {
80 bc_data.push_back(DataFromBc());
81 std::vector<double> mydata;
82 CHKERR bc_data.back().getBcData(mydata, &(*it));
83 CHKERR bc_data.back().getEntitiesFromBc(mField, &(*it));
84 }
85 }
86
88}
89
91
92inline auto get_rotation_from_vector(FTensor1 &t_omega) {
97 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
98 t_R(i, j) = t_kd(i, j);
99
100 const double angle = std::sqrt(t_omega(i) * t_omega(i));
101 if (std::abs(angle) < 1e-18)
102 return t_R;
103
105 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
106 const double a = sin(angle) / angle;
107 const double ss_2 = sin(angle / 2.);
108 const double b = 2. * ss_2 * ss_2 / (angle * angle);
109 t_R(i, j) += a * t_Omega(i, j);
110 t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
111
112 return t_R;
113};
114
115inline auto get_displacement(VectorDouble3 &coords, FTensor1 t_centr,
116 FTensor1 t_normal, double theta) {
120
121 FTensor1 t_omega;
122 FTensor1 t_coords(coords[0], coords[1], coords[2]);
123 const double a = std::sqrt(t_normal(i) * t_normal(i));
124 t_omega(i) = t_normal(i) * (theta / a);
125 auto t_R = get_rotation_from_vector(t_omega);
126 FTensor1 t_delta, t_disp;
127 t_delta(i) = t_centr(i) - t_coords(i);
128 t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
129
130 return VectorDouble({t_disp(0), t_disp(1), t_disp(2)});
131};
132
134 std::vector<DataFromBc> &bc_data) {
135
137
139 if (it->getName().compare(0, 8, "ROTATION") == 0) {
140 bc_data.push_back(DataFromBc());
141 std::vector<double> mydata;
142 CHKERR bc_data.back().getEntitiesFromBc(mField, &(*it));
143 CHKERR it->getAttributes(mydata);
144 if (mydata.size() < 6) {
145 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
146 "6 attributes are required for Rotation (3 center coords + 3 "
147 "angles, (+ 3 optional) flags for xyz)");
148 }
149 for (int ii : {0, 1, 2}) {
150 bc_data.back().bc_flags[ii] = 1;
151 bc_data.back().t_centr(ii) = mydata[ii + 0];
152 bc_data.back().t_normal(ii) = mydata[ii + 3];
153 }
154 if (mydata.size() > 8)
155 for (int ii : {0, 1, 2})
156 bc_data.back().bc_flags[ii] = mydata[6 + ii];
157
158 bc_data.back().scaled_values[0] = bc_data.back().initial_values[0] =
159 bc_data.back().t_normal.l2();
160 bc_data.back().is_rotation = true;
161 }
162 }
163
165}
166
169 DataFromBc &bc_data) {
171 if (bc_data.is_rotation)
172 bc_data.scaled_values = get_displacement(coords, bc_data.t_centr,
173 bc_data.t_normal, bc_data.theta);
175}
176
179 DataFromBc &bc_data) {
181 VectorDouble3 coords(3);
182 CHKERR mField.get_moab().get_coords(&ent, 1, &*coords.begin());
183 CHKERR calculateRotationForDof(coords, bc_data);
185}
186
189
190 bc_data.scaled_values = bc_data.initial_values;
192 bc_data.scaled_values);
193 // for rotation
194 bc_data.theta = bc_data.scaled_values[0];
196}
197
200 if (mapZeroRows.empty() || !methodsOp.empty()) {
201
202 std::vector<DataFromBc> bcData;
205
206 for (auto &bc_it : bcData) {
207
209
210 for (int dim = 0; dim < 3; dim++) {
211
212 auto for_each_dof = [&](auto &dof) {
214 if (dof->getEntType() == MBVERTEX) {
215 CHKERR calculateRotationForDof(dof->getEnt(), bc_it);
216 if (dof->getDofCoeffIdx() == 0 && bc_it.bc_flags[0]) {
217 // set boundary values to field data
218 mapZeroRows[dof->getPetscGlobalDofIdx()] = bc_it.scaled_values[0];
219 dof->getFieldData() = mapZeroRows[dof->getPetscGlobalDofIdx()];
220 }
221 if (dof->getDofCoeffIdx() == 1 && bc_it.bc_flags[1]) {
222 mapZeroRows[dof->getPetscGlobalDofIdx()] = bc_it.scaled_values[1];
223 dof->getFieldData() = mapZeroRows[dof->getPetscGlobalDofIdx()];
224 }
225 if (dof->getDofCoeffIdx() == 2 && bc_it.bc_flags[2]) {
226 mapZeroRows[dof->getPetscGlobalDofIdx()] = bc_it.scaled_values[2];
227 dof->getFieldData() = mapZeroRows[dof->getPetscGlobalDofIdx()];
228 }
229
230 } else {
231 if (dof->getDofCoeffIdx() == 0 && bc_it.bc_flags[0]) {
232 mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
233 dof->getFieldData() = 0;
234 }
235 if (dof->getDofCoeffIdx() == 1 && bc_it.bc_flags[1]) {
236 mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
237 dof->getFieldData() = 0;
238 }
239 if (dof->getDofCoeffIdx() == 2 && bc_it.bc_flags[2]) {
240 mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
241 dof->getFieldData() = 0;
242 }
243 }
245 };
246
249 bc_it.bc_ents[dim], for_each_dof);
250 }
251 }
252 dofsIndices.resize(mapZeroRows.size());
253 dofsValues.resize(mapZeroRows.size());
254 int ii = 0;
255 auto mit = mapZeroRows.begin();
256 for (; mit != mapZeroRows.end(); mit++, ii++) {
257 dofsIndices[ii] = mit->first;
258 dofsValues[ii] = mit->second;
259 }
260 }
262}
263
266
268
269 if (snes_ctx == CTX_SNESNONE) {
270 if (!dofsIndices.empty()) {
272 &*dofsValues.begin(), INSERT_VALUES);
273 }
274 CHKERR VecAssemblyBegin(snes_x);
275 CHKERR VecAssemblyEnd(snes_x);
276 }
277
279}
280
283
284 if (snes_ctx == CTX_SNESNONE) {
285 if (snes_B) {
286 CHKERR MatAssemblyBegin(snes_B, MAT_FINAL_ASSEMBLY);
287 CHKERR MatAssemblyEnd(snes_B, MAT_FINAL_ASSEMBLY);
288 CHKERR MatZeroRowsColumns(snes_B, dofsIndices.size(),
289 dofsIndices.empty() ? PETSC_NULL
290 : &dofsIndices[0],
291 dIag, PETSC_NULL, PETSC_NULL);
292 }
293 if (snes_f) {
294 CHKERR VecAssemblyBegin(snes_f);
295 CHKERR VecAssemblyEnd(snes_f);
296 for (std::vector<int>::iterator vit = dofsIndices.begin();
297 vit != dofsIndices.end(); vit++) {
298 CHKERR VecSetValue(snes_f, *vit, 0, INSERT_VALUES);
299 }
300 CHKERR VecAssemblyBegin(snes_f);
301 CHKERR VecAssemblyEnd(snes_f);
302 }
303 }
304
305 switch (snes_ctx) {
306 case CTX_SNESNONE:
307 break;
308 case CTX_SNESSETFUNCTION: {
309 if (!dofsIndices.empty()) {
310 dofsXValues.resize(dofsIndices.size());
311 const double *a_snes_x;
312 CHKERR VecGetArrayRead(snes_x, &a_snes_x);
313 auto &dofs_by_glob_idx =
315 int idx = 0;
316 for (auto git : dofsIndices) {
317 auto dof_it = dofs_by_glob_idx.find(git);
318 if (dof_it != dofs_by_glob_idx.end()) {
319 dofsXValues[idx] = a_snes_x[dof_it->get()->getPetscLocalDofIdx()];
320 ++idx;
321 } else
322 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
323 "Dof with global %d id not found", git);
324 }
325 CHKERR VecRestoreArrayRead(snes_x, &a_snes_x);
326 }
327 CHKERR VecAssemblyBegin(snes_f);
328 CHKERR VecAssemblyEnd(snes_f);
329
330 if (!dofsIndices.empty()) {
331 int ii = 0;
332 for (std::vector<int>::iterator vit = dofsIndices.begin();
333 vit != dofsIndices.end(); vit++, ii++) {
334 double val = 0;
335 if (!dofsXValues.empty()) {
336 val += dofsXValues[ii];
337 val += -mapZeroRows[*vit]; // in snes it is on the left hand side,
338 // that way -1
339 dofsXValues[ii] = val;
340 }
341 }
343 snes_f, dofsIndices.size(),
344 dofsIndices.empty() ? PETSC_NULL : &*dofsIndices.begin(),
345 dofsXValues.empty() ? PETSC_NULL : &*dofsXValues.begin(),
346 INSERT_VALUES);
347 }
348 CHKERR VecAssemblyBegin(snes_f);
349 CHKERR VecAssemblyEnd(snes_f);
350 } break;
351 case CTX_SNESSETJACOBIAN: {
352
353 CHKERR MatAssemblyBegin(snes_B, MAT_FINAL_ASSEMBLY);
354 CHKERR MatAssemblyEnd(snes_B, MAT_FINAL_ASSEMBLY);
355 CHKERR MatZeroRowsColumns(snes_B, dofsIndices.size(),
356 dofsIndices.empty() ? PETSC_NULL
357 : &*dofsIndices.begin(),
358 dIag, PETSC_NULL, PETSC_NULL);
359
360 } break;
361 default:
362 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown snes stage");
363 }
364
366}
367
370
371 if (mapZeroRows.empty() || !methodsOp.empty()) {
372
373 std::vector<DataFromBc> bcData;
375
376 const FieldEntity_multiIndex *field_ents;
377 CHKERR mField.get_field_ents(&field_ents);
378 auto &field_ents_by_uid = field_ents->get<Unique_mi_tag>();
379
380 VectorDouble3 coords(3);
381
382 for (auto &bc_it : bcData) {
383
385
386 for (int dim = 0; dim < 3; dim++) {
387
388 auto for_each_dof = [&](auto &dof) {
390
391 if (!dim) {
392
393 EntityHandle node = dof->getEnt();
394 if (!dof->getDofCoeffIdx()) {
395 auto eit =
396 field_ents_by_uid.find(FieldEntity::getLocalUniqueIdCalculate(
398 if (eit != field_ents_by_uid.end())
399 noalias(coords) = (*eit)->getEntFieldData();
400 else
401 CHKERR mField.get_moab().get_coords(&node, 1,
402 &*coords.data().begin());
403 }
404 CHKERR calculateRotationForDof(coords, bc_it);
405 if (dof->getDofCoeffIdx() == 0 && bc_it.bc_flags[0]) {
406 mapZeroRows[dof->getPetscGlobalDofIdx()] =
407 coords[0] + bc_it.scaled_values[0];
408 // set boundary values to field data
409 dof->getFieldData() = mapZeroRows[dof->getPetscGlobalDofIdx()];
410 }
411 if (dof->getDofCoeffIdx() == 1 && bc_it.bc_flags[1]) {
412 mapZeroRows[dof->getPetscGlobalDofIdx()] =
413 coords[1] + bc_it.scaled_values[1];
414 dof->getFieldData() = mapZeroRows[dof->getPetscGlobalDofIdx()];
415 }
416 if (dof->getDofCoeffIdx() == 2 && bc_it.bc_flags[2]) {
417 mapZeroRows[dof->getPetscGlobalDofIdx()] =
418 coords[2] + bc_it.scaled_values[2];
419 dof->getFieldData() = mapZeroRows[dof->getPetscGlobalDofIdx()];
420 }
421 } else {
422 if (dof->getDofCoeffIdx() == 0 && bc_it.bc_flags[0]) {
423 mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
424 dof->getFieldData() = 0;
425 }
426 if (dof->getDofCoeffIdx() == 1 && bc_it.bc_flags[1]) {
427 mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
428 dof->getFieldData() = 0;
429 }
430 if (dof->getDofCoeffIdx() == 2 && bc_it.bc_flags[2]) {
431 mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
432 dof->getFieldData() = 0;
433 }
434 }
435
437 };
438
441 bc_it.bc_ents[dim], for_each_dof);
442
443 auto fix_field_dof = [&](auto &dof) {
445 mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
447 };
448
449 for (auto &fix_field : fixFields) {
452 bc_it.bc_ents[dim], fix_field_dof);
453 }
454 }
455 }
456
457 // set vector of values and indices
458 dofsIndices.resize(mapZeroRows.size());
459 dofsValues.resize(mapZeroRows.size());
460 auto mit = mapZeroRows.begin();
461 for (int ii = 0; mit != mapZeroRows.end(); mit++, ii++) {
462 dofsIndices[ii] = mit->first;
463 dofsValues[ii] = mit->second;
464 }
465 }
467}
468
471
472 // function to insert temperature boundary conditions from block
473 auto insert_temp_bc = [&](double &temp, auto &it) {
475 VectorDouble temp_2(1);
476 temp_2[0] = temp;
478 for (int dim = 0; dim < 3; dim++) {
479 Range ents;
480 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), dim, ents,
481 true);
482 if (dim > 1) {
483 Range _edges;
484 CHKERR mField.get_moab().get_adjacencies(ents, 1, false, _edges,
485 moab::Interface::UNION);
486 ents.insert(_edges.begin(), _edges.end());
487 }
488 if (dim > 0) {
489 Range _nodes;
490 CHKERR mField.get_moab().get_connectivity(ents, _nodes, true);
491 ents.insert(_nodes.begin(), _nodes.end());
492 }
493
494 auto for_each_dof = [&](auto &dof) {
496
497 if (dof->getEntType() == MBVERTEX) {
498 mapZeroRows[dof->getPetscGlobalDofIdx()] = temp_2[0];
499 } else {
500 mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
501 }
503 };
504
507 ents, for_each_dof);
508 }
509
511 };
512
513 // Loop over blockset to find the block TEMPERATURE.
514 if (mapZeroRows.empty() || !methodsOp.empty()) {
515 bool is_blockset_defined = false;
516 string blocksetName = "TEMPERATURE";
518 if (it->getName().compare(0, blocksetName.length(), blocksetName) == 0) {
519 std::vector<double> mydata;
520 CHKERR it->getAttributes(mydata);
521
522 if (mydata.empty())
523 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
524 "missing temperature attribute");
525
526 double my_temp = mydata[0];
527 CHKERR insert_temp_bc(my_temp, it);
528
529 is_blockset_defined = true;
530 }
531 }
532
533 // If the block TEMPERATURE is not defined, then
534 // look for temperature boundary conditions
535 if (!is_blockset_defined) {
537 mField, NODESET | TEMPERATURESET, it)) {
539 CHKERR it->getBcDataStructure(mydata);
540 VectorDouble scaled_values(1);
541 scaled_values[0] = mydata.data.value1;
542 CHKERR insert_temp_bc(scaled_values[0], it);
543 }
544 }
545 dofsIndices.resize(mapZeroRows.size());
546 dofsValues.resize(mapZeroRows.size());
547 int ii = 0;
548 std::map<DofIdx, FieldData>::iterator mit = mapZeroRows.begin();
549 for (; mit != mapZeroRows.end(); mit++, ii++) {
550 dofsIndices[ii] = mit->first;
551 dofsValues[ii] = mit->second;
552 }
553 }
554
556}
557
560 if (mapZeroRows.empty()) {
561
562 for (const auto &field_name : fieldNames) {
563
564 auto for_each_dof = [&](auto &dof) {
566 mapZeroRows[dof->getPetscGlobalDofIdx()] = 0;
568 };
569
572 eNts, for_each_dof);
573 }
574
575 dofsIndices.resize(mapZeroRows.size());
576 dofsValues.resize(mapZeroRows.size());
577 int ii = 0;
578 for (auto mit = mapZeroRows.begin(); mit != mapZeroRows.end();
579 ++mit, ++ii) {
580 dofsIndices[ii] = mit->first;
581 dofsValues[ii] = mit->second;
582 }
583 }
585}
586
589
592}
593
596 if (snes_ctx == CTX_SNESNONE) {
597 if (snes_B) {
598 CHKERR MatAssemblyBegin(snes_B, MAT_FINAL_ASSEMBLY);
599 CHKERR MatAssemblyEnd(snes_B, MAT_FINAL_ASSEMBLY);
600 CHKERR MatZeroRowsColumns(snes_B, dofsIndices.size(),
601 dofsIndices.empty() ? PETSC_NULL
602 : &*dofsIndices.begin(),
603 dIag, PETSC_NULL, PETSC_NULL);
604 }
605 if (snes_f) {
606 CHKERR VecAssemblyBegin(snes_f);
607 CHKERR VecAssemblyEnd(snes_f);
608 int ii = 0;
609 for (std::vector<int>::iterator vit = dofsIndices.begin();
610 vit != dofsIndices.end(); vit++, ii++) {
611 CHKERR VecSetValue(snes_f, *vit, dofsValues[ii], INSERT_VALUES);
612 }
613 CHKERR VecAssemblyBegin(snes_f);
614 CHKERR VecAssemblyEnd(snes_f);
615 }
616 }
617
618 switch (snes_ctx) {
619 case CTX_SNESNONE: {
620 } break;
621 case CTX_SNESSETFUNCTION: {
622 CHKERR VecAssemblyBegin(snes_f);
623 CHKERR VecAssemblyEnd(snes_f);
624 int ii = 0;
625 for (std::vector<int>::iterator vit = dofsIndices.begin();
626 vit != dofsIndices.end(); vit++, ii++) {
627 CHKERR VecSetValue(snes_f, *vit, dofsValues[ii], INSERT_VALUES);
628 }
629 CHKERR VecAssemblyBegin(snes_f);
630 CHKERR VecAssemblyEnd(snes_f);
631 } break;
632 case CTX_SNESSETJACOBIAN: {
633 CHKERR MatAssemblyBegin(snes_B, MAT_FINAL_ASSEMBLY);
634 CHKERR MatAssemblyEnd(snes_B, MAT_FINAL_ASSEMBLY);
635 CHKERR MatZeroRowsColumns(snes_B, dofsIndices.size(),
636 dofsIndices.empty() ? PETSC_NULL
637 : &*dofsIndices.begin(),
638 dIag, PETSC_NULL, PETSC_NULL);
639 } break;
640 default:
641 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown snes stage");
642 }
643
645}
646
649
651 auto prb_mng = mField.getInterface<ProblemsManager>();
652 auto field_ptr = mField.get_field_structure(fieldName);
653 const int nb_coefficients = field_ptr->getNbOfCoeffs();
654
655 bcDataPtr = boost::make_shared<vector<DataFromBc>>();
658
659 auto get_dim = [&](const Range &ents) {
660 for (auto d : {3, 2, 1})
661 if (ents.num_of_dimension(d))
662 return d;
663 return 0;
664 };
665
666 auto get_adj_ents = [&](const Range &ents) {
667 Range verts;
668 CHKERR mField.get_moab().get_connectivity(ents, verts, true);
669 const auto dim = get_dim(ents);
670 for (size_t d = 1; d < dim; ++d)
671 CHKERR mField.get_moab().get_adjacencies(ents, d, false, verts,
672 moab::Interface::UNION);
673 verts.merge(ents);
674 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
675 return verts;
676 };
677
678 auto remove_dofs_from_dirichlet_data = [&]() {
680 for (auto &bc_it : *bcDataPtr)
681 for (auto &ents : bc_it.bc_ents)
682 for (int i = 0; i != nb_coefficients; i++)
683 if (bc_it.bc_flags[i])
684 CHKERR prb_mng->removeDofsOnEntities(problemName, fieldName,
685 // get_adj_ents(ents), i, i);
686 ents, i, i);
687
689 };
690
691 auto remove_dofs_from_dirichlet_data_non_distributed = [&]() {
693 for (auto &bc_it : *bcDataPtr)
694 for (auto &ents : bc_it.bc_ents)
695 for (int i = 0; i != nb_coefficients; i++)
696 if (bc_it.bc_flags[i])
697 CHKERR prb_mng->removeDofsOnEntitiesNotDistributed(
698 // problemName, fieldName, get_adj_ents(ents), i, i);
699 problemName, fieldName, ents, i, i, QUIET);
700
702 };
703
704 if (isPartitioned)
705 CHKERR remove_dofs_from_dirichlet_data();
706 else
707 CHKERR remove_dofs_from_dirichlet_data_non_distributed();
708
710}
711
714 if (!bcDataPtr)
716
717 for (auto &bc_it : *bcDataPtr) {
718 Range all_bc_ents;
720 auto method = getEntMethodPtr(bc_it);
721 for (auto &ents : bc_it.bc_ents)
722 all_bc_ents.merge(ents);
723 CHKERR mField.loop_entities(fieldName, *method, &all_bc_ents);
724 }
725
727}
728
730 const MoFEM::CubitMeshSets *it) {
732 // get data structure for boundary condition
733 CHKERR it->getBcDataStructure(mydata);
734 this->scaled_values[0] = mydata.data.value1;
735 this->scaled_values[1] = mydata.data.value2;
736 this->scaled_values[2] = mydata.data.value3;
737 this->initial_values = this->scaled_values;
738 this->bc_flags[0] = (int)mydata.data.flag1;
739 this->bc_flags[1] = (int)mydata.data.flag2;
740 this->bc_flags[2] = (int)mydata.data.flag3;
742}
743
744MoFEMErrorCode DataFromBc::getBcData(std::vector<double> &mydata,
745 const MoFEM::CubitMeshSets *it) {
747 CHKERR it->getAttributes(mydata);
748 if (mydata.size() < 6) {
749 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
750 "six attributes are required for given BC blockset (3 values + "
751 "3 flags)");
752 }
753 for (unsigned int ii = 0; ii < 3; ii++) {
754 this->initial_values[ii] = mydata[ii];
755 this->bc_flags[ii] = (int)mydata[ii + 3];
756 }
759}
760
762 const MoFEM::CubitMeshSets *it) {
764 for (int dim = 0; dim < 3; dim++) {
766 bc_ents[dim], true);
767 if (dim > 1) {
768 Range edges;
769 CHKERR mField.get_moab().get_adjacencies(bc_ents[dim], 1, false, edges,
770 moab::Interface::UNION);
771 bc_ents[1].insert(edges.begin(), edges.end());
772 }
773 if (dim > 0) {
774 Range nodes;
775 CHKERR mField.get_moab().get_connectivity(bc_ents[dim], nodes, true);
776 bc_ents[0].insert(nodes.begin(), nodes.end());
777 }
778 }
780}
781
783
785
786 const Problem *problem_ptr;
787 CHKERR mField.get_problem(problemName.c_str(), &problem_ptr);
788 const double *array;
789 CHKERR VecGetArrayRead(internal, &array);
790
791 auto field_ptr = mField.get_field_structure(fieldName);
792 const int nb_coefficients = field_ptr->getNbOfCoeffs();
793
794 std::vector<int> ghosts(nb_coefficients);
795 for (int g = 0; g != nb_coefficients; ++g)
796 ghosts[g] = g;
797
798 Vec v;
799 CHKERR VecCreateGhost(
800 mField.get_comm(), (mField.get_comm_rank() ? 0 : nb_coefficients),
801 nb_coefficients, (mField.get_comm_rank() ? nb_coefficients : 0),
802 &*ghosts.begin(), &v);
803
806
807 const int id = it->getMeshsetId();
808 VectorDouble &reaction_vec = reactionsMap[id];
809 reaction_vec.resize(nb_coefficients);
810 reaction_vec.clear();
811
812 Range verts;
813 for (int dim = 0; dim != 3; ++dim) {
814 Range ents;
815 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), dim, ents,
816 true);
817 Range nodes;
818 CHKERR mField.get_moab().get_connectivity(ents, nodes, true);
819 verts.insert(nodes.begin(), nodes.end());
820 }
821
822 auto for_each_dof = [&](auto &dof) {
824 reaction_vec[dof->getDofCoeffIdx()] += array[dof->getPetscLocalDofIdx()];
826 };
827
830 verts, for_each_dof);
831
832 double *res_array;
833
834 CHKERR VecGetArray(v, &res_array);
835 for (int dd = 0; dd != reaction_vec.size(); ++dd)
836 res_array[dd] = reaction_vec[dd];
837 CHKERR VecRestoreArray(v, &res_array);
838
839 CHKERR VecGetArray(v, &res_array);
840 for (int dd = 0; dd != reaction_vec.size(); ++dd)
841 reaction_vec[dd] = res_array[dd];
842 CHKERR VecRestoreArray(v, &res_array);
843 }
844
845 CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
846 CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
847 CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
848 CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
849
850 CHKERR VecDestroy(&v);
851 CHKERR VecRestoreArrayRead(internal, &array);
853}
auto get_rotation_from_vector(FTensor1 &t_omega)
Definition: DirichletBC.cpp:92
static MoFEMErrorCode set_numered_dofs_on_ents(const Problem *problem_ptr, const FieldBitNumber bit_number, Range &ents, boost::function< MoFEMErrorCode(const boost::shared_ptr< MoFEM::NumeredDofEntity > &dof)> for_each_dof)
Definition: DirichletBC.cpp:11
auto get_displacement(VectorDouble3 &coords, FTensor1 t_centr, FTensor1 t_normal, double theta)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
constexpr double a
Kronecker Delta class.
@ QUIET
Definition: definitions.h:208
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ TEMPERATURESET
Definition: definitions.h:155
@ NODESET
Definition: definitions.h:146
@ DISPLACEMENTSET
Definition: definitions.h:150
@ BLOCKSET
Definition: definitions.h:148
@ 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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const int dim
@ F
constexpr auto t_kd
void temp(int x, int y=10)
Definition: simple.cpp:4
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, member< FieldEntity, UId, &FieldEntity::localUId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FieldEntity::interface_type_RefEntity, EntityHandle, &FieldEntity::getEnt > > > > FieldEntity_multiIndex
MultiIndex container keeps FieldEntity.
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
virtual MoFEMErrorCode loop_entities(const std::string field_name, EntityMethod &method, Range const *const ents=nullptr, int verb=DEFAULT_VERBOSITY)=0
Loop over field entities.
#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.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
char FieldBitNumber
Field bit number.
Definition: Types.hpp:28
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr auto field_name
constexpr double g
Data from Cubit blocksets.
Definition: DirichletBC.hpp:31
VectorInt bc_flags
Definition: DirichletBC.hpp:34
FTensor::Tensor1< double, 3 > t_centr
Definition: DirichletBC.hpp:40
MoFEMErrorCode getBcData(DisplacementCubitBcData &mydata, const MoFEM::CubitMeshSets *it)
MoFEMErrorCode getEntitiesFromBc(MoFEM::Interface &mField, const MoFEM::CubitMeshSets *it)
VectorDouble scaled_values
Definition: DirichletBC.hpp:32
FTensor::Tensor1< double, 3 > t_normal
Definition: DirichletBC.hpp:39
VectorDouble initial_values
Definition: DirichletBC.hpp:33
double theta
Definition: DirichletBC.hpp:41
bool is_rotation
Definition: DirichletBC.hpp:38
Range bc_ents[3]
Definition: DirichletBC.hpp:35
DirichletDisplacementBc(MoFEM::Interface &m_field, const std::string &field_name, Mat Aij, Vec X, Vec F, string blockset_name="DISPLACEMENT")
Definition: DirichletBC.cpp:37
const std::string blocksetName
Definition: DirichletBC.hpp:74
MoFEMErrorCode postProcess()
function is run at the end of loop
std::map< DofIdx, FieldData > mapZeroRows
Definition: DirichletBC.hpp:70
MoFEMErrorCode getRotationBcFromBlock(std::vector< DataFromBc > &bc_data)
Get the Rotation Bc From Block object Use ROTATION blockset name with 7 atributes: 1,...
const std::string fieldName
field name to set Dirichlet BC
Definition: DirichletBC.hpp:60
double dIag
diagonal value set on zeroed column and rows
Definition: DirichletBC.hpp:61
std::vector< int > dofsIndices
Definition: DirichletBC.hpp:71
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode getBcDataFromSetsAndBlocks(std::vector< DataFromBc > &bc_data)
Get the Bc Data From Sets And Blocks object Use DISPLACEMENT blockset name (default) with 6 atributes...
Definition: DirichletBC.cpp:64
MoFEMErrorCode applyScaleBcData(DataFromBc &bc_data)
std::vector< double > dofsXValues
Definition: DirichletBC.hpp:73
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: DirichletBC.hpp:76
std::vector< double > dofsValues
Definition: DirichletBC.hpp:72
MoFEM::Interface & mField
Definition: DirichletBC.hpp:59
virtual MoFEMErrorCode iNitialize()
MoFEMErrorCode calculateRotationForDof(VectorDouble3 &coords, DataFromBc &bc_data)
Calculate displacements from rotation for particular dof.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
virtual boost::shared_ptr< EntityMethod > getEntMethodPtr(DataFromBc &data)
boost::shared_ptr< vector< DataFromBc > > bcDataPtr
std::vector< std::string > fieldNames
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode iNitialize()
std::vector< std::string > fixFields
MoFEMErrorCode iNitialize()
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
const Problem * problemPtr
raw pointer to problem
Managing BitRefLevels.
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
this struct keeps basic methods for moab meshset about material and boundary conditions
MoFEMErrorCode getAttributes(std::vector< double > &attributes) const
get Cubit block attributes
MoFEMErrorCode getBcDataStructure(CUBIT_BC_DATA_TYPE &data) const
MoFEMErrorCode getMeshsetIdEntitiesByDimension(Interface &moab, const int dimension, Range &entities, const bool recursive=false) const
get entities form meshset
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:72
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
FieldCoefficientsNumber getNbOfCoeffs() const
Get number of field coefficients.
keeps basic data about problem
auto & getNumeredColDofsPtr() const
get access to numeredColDofsPtr storing DOFs on cols
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
Vec & snes_f
residual
Vec & snes_x
state vector
Mat & snes_B
preconditioner of jacobian matrix
SNESContext snes_ctx
Vec & ts_F
residual vector
Vec & ts_u
state vector
Mat & ts_B
Preconditioner for ts_A.
Definition of the temperature bc data structure.
Definition: BCData.hpp:301
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
std::string problemName
MoFEM::Interface & mField
ReactionsMap reactionsMap
MoFEMErrorCode calculateReactions(Vec &internal)
calculate reactions from a given vector
std::string fieldName