8#ifndef _MIX_TRANPORT_ELEMENT_HPP_
9#define _MIX_TRANPORT_ELEMENT_HPP_
93 ids.empty() ? PETSC_NULL : &ids[0],
94 PETSC_COPY_VALUES, &is_local);
95 CHKERR ISAllGather(is_local, is);
96 CHKERR ISDestroy(&is_local);
110 const double y,
const double z,
127 const double y,
const double z,
128 MatrixDouble3by3 &inv_k) {
131 for (
int dd = 0; dd < 3; dd++) {
147 const double x,
const double y,
148 const double z,
double &value) {
164 const double y,
const double z,
179 std::map<int, BlockData>
189 MoFEMErrorCode
addFields(
const std::string &values,
const std::string &fluxes,
209 const std::string &values_name) {
244 Mat_Thermal temp_data;
245 CHKERR it->getAttributeDataStructure(temp_data);
247 temp_data.data.Conductivity;
248 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
250 it->meshset, MBTET,
setOfBlocks[it->getMeshsetId()].tEts,
true);
252 setOfBlocks[it->getMeshsetId()].tEts, MBTET,
"MIX");
256 setOfBlocks[it->getMeshsetId()].tEts, 2,
false, skeleton,
257 moab::Interface::UNION);
299 BitRefLevel().set(0), BitRefLevel().set(), MBTET, done_tets);
306 ref_level, BitRefLevel().set(), MBTET, ref_tets);
307 ref_tets = subtract(ref_tets, done_tets);
313 BitRefLevel().set(0), BitRefLevel().set(), MBTRI, done_faces);
320 ref_level, BitRefLevel().set(), MBTRI, ref_faces);
321 ref_faces = subtract(ref_faces, done_faces);
340 ProblemsManager *prb_mng_ptr;
342 CHKERR prb_mng_ptr->buildProblem(
"MIX",
true);
345 CHKERR prb_mng_ptr->partitionProblem(
"MIX");
346 CHKERR prb_mng_ptr->partitionFiniteElements(
"MIX");
348 CHKERR prb_mng_ptr->partitionGhostDofs(
"MIX");
357 std::vector<EntityHandle> &map_gauss_pts)
362 EntitiesFieldData::EntData &data) {
370 double *error_flux_ptr;
372 th_error_flux, &fe_ent, 1, (
const void **)&error_flux_ptr);
377 double *error_div_ptr;
379 th_error_div, &fe_ent, 1, (
const void **)&error_div_ptr);
384 double *error_jump_ptr;
386 th_error_jump, &fe_ent, 1, (
const void **)&error_jump_ptr);
392 "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
393 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
394 for (vector<EntityHandle>::iterator vit =
mapGaussPts.begin();
402 "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
403 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
404 for (vector<EntityHandle>::iterator vit =
mapGaussPts.begin();
412 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
413 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
414 for (vector<EntityHandle>::iterator vit =
mapGaussPts.begin();
431 PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
434 CHKERR AddHOOps<3, 3, 3>::add(post_proc.getOpPtrVector(), {HDIV, L2});
436 auto values_ptr = boost::make_shared<VectorDouble>();
437 auto grad_ptr = boost::make_shared<MatrixDouble>();
438 auto flux_ptr = boost::make_shared<MatrixDouble>();
440 post_proc.getOpPtrVector().push_back(
441 new OpCalculateScalarFieldValues(
"VALUES", values_ptr));
442 post_proc.getOpPtrVector().push_back(
443 new OpCalculateScalarFieldGradient<3>(
"VALUES", grad_ptr));
444 post_proc.getOpPtrVector().push_back(
445 new OpCalculateHVecVectorField<3>(
"FLUXES", flux_ptr));
447 using OpPPMap = OpPostProcMapInMoab<3, 3>;
449 post_proc.getOpPtrVector().push_back(
451 new OpPPMap(post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
453 {{
"VALUES", values_ptr}},
455 {{
"GRAD", grad_ptr}, {
"FLUXES", flux_ptr}},
461 post_proc.getOpPtrVector().push_back(
new OpPostProc(
462 post_proc.getPostProcMesh(), post_proc.getMapGaussPts()));
466 CHKERR post_proc.writeFile(out_file.c_str());
477 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(
"MIX", &Aij);
490 CHKERR MatZeroEntries(Aij);
492 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
493 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
494 CHKERR VecZeroEntries(D0);
495 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
496 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
498 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
499 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
502 "MIX",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
515 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
516 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
517 CHKERR VecAssemblyBegin(D0);
518 CHKERR VecAssemblyEnd(D0);
542 CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
543 CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
544 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
545 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
551 CHKERR VecNorm(
F, NORM_2, &nrm2_F);
552 PetscPrintf(PETSC_COMM_WORLD,
"nrm2_F = %6.4e\n", nrm2_F);
562 CHKERR getDirichletBCIndices(&essential_bc_ids);
563 CHKERR MatZeroRowsColumnsIS(Aij, essential_bc_ids, 1, D0,
F);
564 CHKERR ISDestroy(&essential_bc_ids);
574 CHKERR VecNorm(
F, NORM_2, &nrm2_F);
575 PetscPrintf(PETSC_COMM_WORLD,
"With BC nrm2_F = %6.4e\n", nrm2_F);
588 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
589 CHKERR KSPSetOperators(solver, Aij, Aij);
590 CHKERR KSPSetFromOptions(solver);
593 CHKERR KSPDestroy(&solver);
597 CHKERR VecNorm(
D, NORM_2, &nrm2_D);
599 PetscPrintf(PETSC_COMM_WORLD,
"nrm2_D = %6.4e\n", nrm2_D);
601 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
602 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
606 "MIX",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
615 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
616 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
636 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
637 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
643 std::vector<int> ids;
644 ids.insert(ids.begin(), bcIndices.begin(), bcIndices.end());
645 std::vector<double> vals(ids.size(), 0);
646 CHKERR VecSetValues(
F, ids.size(), &*ids.begin(), &*vals.begin(),
653 CHKERR VecNorm(
F, NORM_2, &nrm2_F);
654 PetscPrintf(PETSC_COMM_WORLD,
"nrm2_F = %6.4e\n", nrm2_F);
655 const double eps = 1e-8;
689 const Problem *problem_ptr;
692 "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
693 "jump^2 = %6.4e error tot^2 = %6.4e\n",
694 problem_ptr->getNbDofsRow(), sumErrorFlux, sumErrorDiv,
695 sumErrorJump, sumErrorFlux + sumErrorDiv + sumErrorJump);
727 const std::string flux_name, Mat aij, Vec f)
729 flux_name, flux_name,
731 cTx(ctx), Aij(aij),
F(f) {
735 MatrixDouble
NN, transNN;
751 EntitiesFieldData::EntData &row_data,
752 EntitiesFieldData::EntData &col_data) {
754 if (Aij == PETSC_NULL)
756 if (row_data.getIndices().size() == 0)
758 if (col_data.getIndices().size() == 0)
760 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
761 int nb_row = row_data.getIndices().size();
762 int nb_col = col_data.getIndices().size();
763 NN.resize(nb_row, nb_col,
false);
767 invK.resize(3, 3,
false);
770 &invK(0, 0), &invK(0, 1), &invK(0, 2), &invK(1, 0), &invK(1, 1),
771 &invK(1, 2), &invK(2, 0), &invK(2, 1), &invK(2, 2));
773 auto t_n_hdiv_row = row_data.getFTensor1N<3>();
775 int nb_gauss_pts = row_data.getN().size1();
776 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
778 double w = getGaussPts()(3, gg) * getVolume();
779 const double x = getCoordsAtGaussPts()(gg, 0);
780 const double y = getCoordsAtGaussPts()(gg, 1);
781 const double z = getCoordsAtGaussPts()(gg, 2);
784 for (
int kk = 0; kk != nb_row; kk++) {
786 &col_data.getVectorN<3>(gg)(0,
HVEC0),
787 &col_data.getVectorN<3>(gg)(0,
HVEC1),
788 &col_data.getVectorN<3>(gg)(0,
HVEC2), 3);
789 t_row(
j) = w * t_n_hdiv_row(
i) * t_inv_k(
i,
j);
790 for (
int ll = 0; ll != nb_col; ll++) {
791 NN(kk, ll) += t_row(
j) * t_n_hdiv_col(
j);
797 Mat
a = (Aij != PETSC_NULL) ? Aij : getFEMethod()->ts_B;
798 CHKERR MatSetValues(
a, nb_row, &row_data.getIndices()[0], nb_col,
799 &col_data.getIndices()[0], &NN(0, 0), ADD_VALUES);
801 if (row_side != col_side || row_type != col_type) {
802 transNN.resize(nb_col, nb_row);
803 noalias(transNN) = trans(NN);
804 CHKERR MatSetValues(
a, nb_col, &col_data.getIndices()[0], nb_row,
805 &row_data.getIndices()[0], &transNN(0, 0),
820 EntitiesFieldData::EntData &data) {
824 int nb_row = data.getIndices().size();
828 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
833 invK.resize(3, 3,
false);
838 &invK(0, 0), &invK(0, 1), &invK(0, 2), &invK(1, 0), &invK(1, 1),
839 &invK(1, 2), &invK(2, 0), &invK(2, 1), &invK(2, 2));
841 auto t_n_hdiv = data.getFTensor1N<3>();
843 int nb_gauss_pts = data.getN().size1();
844 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
845 double w = getGaussPts()(3, gg) * getVolume();
846 const double x = getCoordsAtGaussPts()(gg, 0);
847 const double y = getCoordsAtGaussPts()(gg, 1);
848 const double z = getCoordsAtGaussPts()(gg, 2);
850 for (
int ll = 0; ll != nb_row; ll++) {
851 Nf[ll] += w * t_n_hdiv(
i) * t_inv_k(
i,
j) * t_flux(
j);
856 CHKERR VecSetValues(
F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
872 string val_name_col, Vec f)
885 EntitiesFieldData::EntData &data) {
887 if (data.getFieldData().size() == 0)
889 int nb_row = data.getIndices().size();
892 divVec.resize(data.getN().size2() / 3, 0);
893 if (divVec.size() != data.getIndices().size()) {
895 "data inconsistency");
897 int nb_gauss_pts = data.getN().size1();
900 auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
903 for (; gg < nb_gauss_pts; gg++) {
904 double w = getGaussPts()(3, gg) * getVolume();
905 for (
auto &
v : divVec) {
906 v = t_base_diff_hdiv(
i,
i);
911 CHKERR VecSetValues(
F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
932 string flux_name_col, Mat aij, Vec f)
934 val_name_row, flux_name_col,
936 cTx(ctx), Aij(aij),
F(f) {
944 MatrixDouble
NN, transNN;
961 EntitiesFieldData::EntData &row_data,
962 EntitiesFieldData::EntData &col_data) {
964 if (Aij == PETSC_NULL)
966 if (row_data.getFieldData().size() == 0)
968 if (col_data.getFieldData().size() == 0)
970 int nb_row = row_data.getFieldData().size();
971 int nb_col = col_data.getFieldData().size();
972 NN.resize(nb_row, nb_col);
974 divVec.resize(nb_col,
false);
977 auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
979 int nb_gauss_pts = row_data.getN().size1();
980 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
981 double w = getGaussPts()(3, gg) * getVolume();
982 for (
auto &
v : divVec) {
983 v = t_base_diff_hdiv(
i,
i);
986 noalias(NN) += w * outer_prod(row_data.getN(gg), divVec);
988 CHKERR MatSetValues(Aij, nb_row, &row_data.getIndices()[0], nb_col,
989 &col_data.getIndices()[0], &NN(0, 0), ADD_VALUES);
990 transNN.resize(nb_col, nb_row);
991 ublas::noalias(transNN) = -trans(NN);
992 CHKERR MatSetValues(Aij, nb_col, &col_data.getIndices()[0], nb_row,
993 &row_data.getIndices()[0], &transNN(0, 0),
999 EntitiesFieldData::EntData &data) {
1001 if (data.getIndices().size() == 0)
1003 if (data.getIndices().size() != data.getN().size2()) {
1005 "data inconsistency");
1007 int nb_row = data.getIndices().size();
1010 int nb_gauss_pts = data.getN().size1();
1011 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1012 double w = getGaussPts()(3, gg) * getVolume();
1015 CHKERR VecSetValues(
F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
1036 EntitiesFieldData::EntData &data) {
1038 if (data.getFieldData().size() == 0)
1040 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1041 int nb_row = data.getFieldData().size();
1042 Nf.resize(nb_row,
false);
1044 int nb_gauss_pts = data.getN().size1();
1045 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1046 double w = getGaussPts()(3, gg) * getVolume();
1047 const double x = getCoordsAtGaussPts()(gg, 0);
1048 const double y = getCoordsAtGaussPts()(gg, 1);
1049 const double z = getCoordsAtGaussPts()(gg, 2);
1052 noalias(Nf) += w * data.getN(gg) * flux;
1054 CHKERR VecSetValues(
F, nb_row, &data.getIndices()[0], &Nf[0], ADD_VALUES);
1093 EntitiesFieldData::EntData &data) {
1095 if (data.getFieldData().size() == 0)
1097 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1098 nF.resize(data.getIndices().size());
1100 int nb_gauss_pts = data.getN().size1();
1101 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1102 const double x = getCoordsAtGaussPts()(gg, 0);
1103 const double y = getCoordsAtGaussPts()(gg, 1);
1104 const double z = getCoordsAtGaussPts()(gg, 2);
1109 double w = getGaussPts()(2, gg) * 0.5;
1110 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1112 w * prod(data.getVectorN<3>(gg), getNormalsAtGaussPts(gg)) *
1115 noalias(nF) += w * prod(data.getVectorN<3>(gg), getNormal()) * value;
1118 Vec f = (
F != PETSC_NULL) ?
F : getFEMethod()->ts_F;
1119 CHKERR VecSetValues(f, data.getIndices().size(), &data.getIndices()[0],
1120 &nF[0], ADD_VALUES);
1152 EntitiesFieldData::EntData &data) {
1154 if (data.getFieldData().size() == 0)
1156 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1157 int nb_dofs = data.getFieldData().size();
1158 int nb_gauss_pts = data.getN().size1();
1159 if (3 * nb_dofs !=
static_cast<int>(data.getN().size2())) {
1161 "wrong number of dofs");
1163 NN.resize(nb_dofs, nb_dofs);
1172 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1174 normal_ptr = &getNormalsAtGaussPts(0)[0];
1177 normal_ptr = &getNormal()[0];
1183 auto t_n_hdiv_row = data.getFTensor1N<3>();
1188 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1191 const double x = getCoordsAtGaussPts()(gg, 0);
1192 const double y = getCoordsAtGaussPts()(gg, 1);
1193 const double z = getCoordsAtGaussPts()(gg, 2);
1200 double w = getGaussPts()(2, gg);
1202 nrm2 = sqrt(t_normal(
i) * t_normal(
i));
1207 for (
int ll = 0; ll != nb_dofs; ll++) {
1210 &data.getVectorN<3>(gg)(0,
HVEC0),
1211 &data.getVectorN<3>(gg)(0,
HVEC1),
1212 &data.getVectorN<3>(gg)(0,
HVEC2), 3);
1213 for (
int kk = 0; kk <= ll; kk++) {
1214 NN(ll, kk) += w * t_n_hdiv_row(
i) * t_n_hdiv_col(
i);
1218 Nf[ll] += w * t_n_hdiv_row(
i) * t_normal(
i) * flux / nrm2;
1223 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1225 nrm2 = sqrt(t_normal(
i) * t_normal(
i));
1229 cTx.
bcIndices.insert(data.getIndices().begin(), data.getIndices().end());
1236 CHKERR VecSetOption(cTx.
D0, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1237 CHKERR VecSetValues(cTx.
D0, nb_dofs, &*data.getIndices().begin(),
1238 &*Nf.begin(), INSERT_VALUES);
1239 for (
int dd = 0; dd != nb_dofs; ++dd)
1240 data.getFieldDofs()[dd]->getFieldData() = Nf[dd];
1255 const std::string val_name =
"VALUES")
1263 EntitiesFieldData::EntData &data) {
1265 if (data.getFieldData().size() == 0)
1268 int nb_gauss_pts = data.getN().size1();
1270 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1272 inner_prod(trans(data.getN(gg)), data.getFieldData());
1288 const std::string val_name =
"VALUES")
1295 EntitiesFieldData::EntData &data) {
1297 if (data.getFieldData().size() == 0)
1299 int nb_gauss_pts = data.getDiffN().size1();
1301 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1302 ublas::matrix_column<MatrixDouble> values_grad_at_gauss_pts(
1304 noalias(values_grad_at_gauss_pts) =
1305 prod(trans(data.getDiffN(gg)), data.getFieldData());
1327 EntitiesFieldData::EntData &data) {
1330 if (data.getFieldData().size() == 0)
1332 int nb_gauss_pts = data.getDiffN().size1();
1333 int nb_dofs = data.getFieldData().size();
1336 if (type == MBTRI && side == 0) {
1340 divVec.resize(nb_dofs);
1343 auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
1345 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1346 for (
auto &
v : divVec) {
1347 v = t_base_diff_hdiv(
i,
i);
1351 ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1354 prod(trans(data.getVectorN<3>(gg)), data.getFieldData());
1382 EntitiesFieldData::EntData &data) {
1386 invK.resize(3, 3,
false);
1387 int nb_gauss_pts = data.getN().size1();
1388 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1390 Tag th_error_flux, th_error_div;
1392 "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
1393 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1394 double *error_flux_ptr;
1396 th_error_flux, &fe_ent, 1, (
const void **)&error_flux_ptr);
1399 "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
1400 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1401 double *error_div_ptr;
1403 th_error_div, &fe_ent, 1, (
const void **)&error_div_ptr);
1407 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1408 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1409 double *error_jump_ptr;
1411 th_error_jump, &fe_ent, 1, (
const void **)&error_jump_ptr);
1412 *error_jump_ptr = 0;
1415 const double h = pow(getVolume() * 12 / std::sqrt(2), (
double)1 / 3);
1417 for (
int ff = 0; ff != 4; ff++) {
1420 double *error_face_jump_ptr;
1422 th_error_jump, &face, 1, (
const void **)&error_face_jump_ptr);
1423 *error_face_jump_ptr = (1 / sqrt(
h)) * sqrt(*error_face_jump_ptr);
1424 *error_face_jump_ptr = pow(*error_face_jump_ptr, 2);
1425 *error_jump_ptr += *error_face_jump_ptr;
1428 *error_flux_ptr = 0;
1430 deltaFlux.resize(3,
false);
1431 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1432 double w = getGaussPts()(3, gg) * getVolume();
1433 const double x = getCoordsAtGaussPts()(gg, 0);
1434 const double y = getCoordsAtGaussPts()(gg, 1);
1435 const double z = getCoordsAtGaussPts()(gg, 2);
1442 ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1444 ublas::matrix_column<MatrixDouble> values_grad_at_gauss_pts(
1446 noalias(deltaFlux) =
1447 prod(invK, flux_at_gauss_pt) + values_grad_at_gauss_pts;
1448 *error_flux_ptr += w * inner_prod(deltaFlux, deltaFlux);
1450 *error_div_ptr =
h * sqrt(*error_div_ptr);
1451 *error_div_ptr = pow(*error_div_ptr, 2);
1452 cTx.
errorMap[sqrt(*error_flux_ptr + *error_div_ptr + *error_jump_ptr)] =
1459 *error_jump_ptr * getVolume();
1484 :
public VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator {
1491 EntitiesFieldData::EntData &data) {
1493 if (data.getFieldData().size() == 0)
1495 int nb_gauss_pts = data.getN().size1();
1496 valMap[getFaceSense()].resize(nb_gauss_pts);
1497 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1498 valMap[getFaceSense()][gg] =
1499 inner_prod(trans(data.getN(gg)), data.getFieldData());
1510 volSideFe(ctx.mField), cTx(ctx) {
1515 EntitiesFieldData::EntData &data) {
1518 if (type == MBTRI) {
1519 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1524 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1525 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1526 double *error_jump_ptr;
1528 th_error_jump, &fe_ent, 1, (
const void **)&error_jump_ptr);
1529 *error_jump_ptr = 0;
1534 if (cTx.
mField.
get_moab().contains_entities(essential_bc_meshset,
1542 CHKERR loopSideVolumes(
"MIX", volSideFe);
1545 int nb_gauss_pts = data.getN().size1();
1548 if (valMap.size() == 1) {
1549 if (
static_cast<int>(valMap.begin()->second.size()) != nb_gauss_pts) {
1551 "wrong number of integration points");
1553 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
1554 const double x = getCoordsAtGaussPts()(gg, 0);
1555 const double y = getCoordsAtGaussPts()(gg, 1);
1556 const double z = getCoordsAtGaussPts()(gg, 2);
1559 double w = getGaussPts()(2, gg);
1560 if (
static_cast<int>(getNormalsAtGaussPts().size1()) ==
1562 w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1566 *error_jump_ptr += w * pow(value - valMap.begin()->second[gg], 2);
1568 }
else if (valMap.size() == 2) {
1569 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
1570 double w = getGaussPts()(2, gg);
1571 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1572 w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1576 double delta = valMap.at(1)[gg] - valMap.at(-1)[gg];
1577 *error_jump_ptr += w * pow(
delta, 2);
1581 "data inconsistency, wrong number of neighbors "
1582 "valMap.size() = %d",
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define BITREFLEVEL_SIZE
max number of refinements
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MAT_THERMALSET
block name is "MAT_THERMAL"
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
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_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
implementation of Data Operators for Forces and Sources
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
static constexpr double delta
data for evaluation of het conductivity and heat capacity elements
Range tEts
constatins elements in block set
definition of surface element
MyTriFE(MoFEM::Interface &m_field)
definition of volume element
MyVolumeFE(MoFEM::Interface &m_field)
MixTransportElement & cTx
OpDivTauU_HdivL2(MixTransportElement &ctx, const std::string flux_name_row, string val_name_col, Vec f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate error evaluator
OpError(MixTransportElement &ctx, const std::string field_name)
MixTransportElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Evaluate boundary conditions on fluxes.
FTensor::Index< 'i', 3 > i
MixTransportElement & cTx
OpEvaluateBcOnFluxes(MixTransportElement &ctx, const std::string flux_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate flux at integration point
OpFluxDivergenceAtGaussPts(MixTransportElement &ctx, const std::string field_name)
MixTransportElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate source therms, i.e. .
OpL2Source(MixTransportElement &ctx, const std::string val_name, Vec f)
MixTransportElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpPostProc(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts)
MixTransportElement & cTx
OpRhsBcOnValues(MixTransportElement &ctx, const std::string fluxes_name, Vec f)
Constructor.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
calculate values on adjacent tetrahedra to face
map< int, VectorDouble > & valMap
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpVolSide(map< int, VectorDouble > &val_map)
calculate jump on entities
OpSkeleton(MixTransportElement &ctx, const std::string flux_name)
map< int, VectorDouble > valMap
VolumeElementForcesAndSourcesCoreOnSide volSideFe
volume element to get values from adjacent tets to face
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MixTransportElement & cTx
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
MixTransportElement & cTx
OpTauDotSigma_HdivHdiv(MixTransportElement &ctx, const std::string flux_name, Mat aij, Vec f)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Assemble matrix.
MixTransportElement & cTx
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
OpVDivSigma_L2Hdiv(MixTransportElement &ctx, const std::string val_name_row, string flux_name_col, Mat aij, Vec f)
Constructor.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate values at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpValuesAtGaussPts(MixTransportElement &ctx, const std::string val_name="VALUES")
MixTransportElement & cTx
virtual ~OpValuesAtGaussPts()
Calculate gradients of values at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpValuesGradientAtGaussPts(MixTransportElement &ctx, const std::string val_name="VALUES")
MixTransportElement & cTx
virtual ~OpValuesGradientAtGaussPts()
MoFEMErrorCode postProc(const string out_file)
Post process results.
MoFEMErrorCode destroyMatrices()
destroy matrices
MoFEMErrorCode evaluateError()
Calculate error on elements.
MoFEMErrorCode solveLinearProblem()
solve problem
MoFEMErrorCode getDirichletBCIndices(IS *is)
get dof indices where essential boundary conditions are applied
MoFEMErrorCode createMatrices()
create matrices
MoFEM::Interface & mField
MatrixDouble valuesGradientAtGaussPts
gradients at integration points on element
virtual MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
MoFEMErrorCode calculateResidual()
calculate residual
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
VectorDouble valuesAtGaussPts
values at integration points on element
map< double, EntityHandle > errorMap
MyTriFE feTri
Instance of surface element.
virtual MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
evaluate natural (Dirichlet) boundary conditions
MyVolumeFE feVol
Instance of volume element.
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
Add fields to database.
virtual MoFEMErrorCode getResistivity(const EntityHandle ent, const double x, const double y, const double z, MatrixDouble3by3 &inv_k)
natural (Dirichlet) boundary conditions (set values)
MixTransportElement(MoFEM::Interface &m_field)
construction of this data structure
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
VectorDouble divergenceAtGaussPts
divergence at integration points on element
virtual ~MixTransportElement()
destructor
virtual MoFEMErrorCode getSource(const EntityHandle ent, const double x, const double y, const double z, double &flux)
set source term
MatrixDouble fluxesAtGaussPts
fluxes at integration points on element
MoFEMErrorCode buildProblem(BitRefLevel &ref_level)
Build problem.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
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.
default operator for TRI element
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
@ OPCOL
operator doWork function is executed on FE columns
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Post post-proc data at points from hash maps.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCore * getVolumeFE() const
return pointer to Generic Volume Finite Element object
Volume finite element base.
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)