23 #ifndef _MIX_TRANPORT_ELEMENT_HPP_
24 #define _MIX_TRANPORT_ELEMENT_HPP_
104 std::vector<int> ids;
108 ids.empty() ? PETSC_NULL : &ids[0],
109 PETSC_COPY_VALUES, &is_local);
110 CHKERR ISAllGather(is_local, is);
111 CHKERR ISDestroy(&is_local);
125 const double y,
const double z,
142 const double y,
const double z,
146 for (
int dd = 0;
dd < 3;
dd++) {
162 const double x,
const double y,
163 const double z,
double &value) {
179 const double y,
const double z,
194 std::map<int, BlockData>
224 const std::string &fluxes_name,
const std::string &values_name,
225 const std::string mesh_nodals_positions =
"MESH_NODE_POSITIONS") {
262 mesh_nodals_positions);
264 mesh_nodals_positions);
274 Mat_Thermal temp_data;
275 CHKERR it->getAttributeDataStructure(temp_data);
277 temp_data.data.Conductivity;
278 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
280 it->meshset, MBTET,
setOfBlocks[it->getMeshsetId()].tEts,
true);
282 setOfBlocks[it->getMeshsetId()].tEts, MBTET,
"MIX");
286 setOfBlocks[it->getMeshsetId()].tEts, 2,
false, skeleton,
287 moab::Interface::UNION);
304 mesh_nodals_positions);
319 mesh_nodals_positions);
346 ref_tets = subtract(ref_tets, done_tets);
360 ref_faces = subtract(ref_faces, done_faces);
379 ProblemsManager *prb_mng_ptr;
381 CHKERR prb_mng_ptr->buildProblem(
"MIX",
true);
384 CHKERR prb_mng_ptr->partitionProblem(
"MIX");
385 CHKERR prb_mng_ptr->partitionFiniteElements(
"MIX");
387 CHKERR prb_mng_ptr->partitionGhostDofs(
"MIX");
396 std::vector<EntityHandle> &map_gauss_pts)
409 double *error_flux_ptr;
411 th_error_flux, &fe_ent, 1, (
const void **)&error_flux_ptr);
416 double *error_div_ptr;
418 th_error_div, &fe_ent, 1, (
const void **)&error_div_ptr);
423 double *error_jump_ptr;
425 th_error_jump, &fe_ent, 1, (
const void **)&error_jump_ptr);
431 "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
432 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
433 for (vector<EntityHandle>::iterator vit =
mapGaussPts.begin();
441 "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
442 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
443 for (vector<EntityHandle>::iterator vit =
mapGaussPts.begin();
451 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
452 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
453 for (vector<EntityHandle>::iterator vit =
mapGaussPts.begin();
489 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(
"MIX", &
Aij);
504 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
505 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
507 CHKERR VecGhostUpdateBegin(
D0, INSERT_VALUES, SCATTER_FORWARD);
508 CHKERR VecGhostUpdateEnd(
D0, INSERT_VALUES, SCATTER_FORWARD);
510 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
511 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
514 "MIX",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
522 feTri.getOpPtrVector().clear();
526 CHKERR VecGhostUpdateBegin(
D0, INSERT_VALUES, SCATTER_REVERSE);
527 CHKERR VecGhostUpdateEnd(
D0, INSERT_VALUES, SCATTER_REVERSE);
547 feTri.getOpPtrVector().clear();
553 CHKERR MatAssemblyBegin(
Aij, MAT_FINAL_ASSEMBLY);
554 CHKERR MatAssemblyEnd(
Aij, MAT_FINAL_ASSEMBLY);
555 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
556 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
562 CHKERR VecNorm(
F, NORM_2, &nrm2_F);
563 PetscPrintf(PETSC_COMM_WORLD,
"nrm2_F = %6.4e\n", nrm2_F);
574 CHKERR MatZeroRowsColumnsIS(
Aij, essential_bc_ids, 1,
D0,
F);
575 CHKERR ISDestroy(&essential_bc_ids);
585 CHKERR VecNorm(
F, NORM_2, &nrm2_F);
586 PetscPrintf(PETSC_COMM_WORLD,
"With BC nrm2_F = %6.4e\n", nrm2_F);
599 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
601 CHKERR KSPSetFromOptions(solver);
604 CHKERR KSPDestroy(&solver);
608 CHKERR VecNorm(
D, NORM_2, &nrm2_D);
610 PetscPrintf(PETSC_COMM_WORLD,
"nrm2_D = %6.4e\n", nrm2_D);
612 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
613 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
617 "MIX",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
626 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
627 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
643 feTri.getOpPtrVector().clear();
646 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
647 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
653 std::vector<int> ids;
655 std::vector<double> vals(ids.size(), 0);
663 CHKERR VecNorm(
F, NORM_2, &nrm2_F);
664 PetscPrintf(PETSC_COMM_WORLD,
"nrm2_F = %6.4e\n", nrm2_F);
665 const double eps = 1e-8;
686 feTri.getOpPtrVector().clear();
698 const Problem *problem_ptr;
701 "Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error "
702 "jump^2 = %6.4e error tot^2 = %6.4e\n",
736 const std::string flux_name, Mat aij,
Vec f)
738 flux_name, flux_name,
763 if (
Aij == PETSC_NULL)
765 if (row_data.getIndices().size() == 0)
767 if (col_data.getIndices().size() == 0)
770 int nb_row = row_data.getIndices().size();
771 int nb_col = col_data.getIndices().size();
772 NN.resize(nb_row, nb_col,
false);
776 invK.resize(3, 3,
false);
779 &
invK(0, 0), &
invK(0, 1), &
invK(0, 2), &
invK(1, 0), &
invK(1, 1),
782 auto t_n_hdiv_row = row_data.getFTensor1N<3>();
784 int nb_gauss_pts = row_data.getN().size1();
785 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
798 for (
int kk = 0; kk != nb_row; kk++) {
800 &col_data.getVectorN<3>(gg)(0,
HVEC0),
801 &col_data.getVectorN<3>(gg)(0,
HVEC1),
802 &col_data.getVectorN<3>(gg)(0,
HVEC2), 3);
803 t_row(
j) =
w * t_n_hdiv_row(
i) * t_inv_k(
i,
j);
804 for (
int ll = 0; ll != nb_col; ll++) {
805 NN(kk, ll) += t_row(
j) * t_n_hdiv_col(
j);
813 &col_data.getIndices()[0], &
NN(0, 0), ADD_VALUES);
815 if (row_side != col_side || row_type != col_type) {
816 transNN.resize(nb_col, nb_row);
819 &row_data.getIndices()[0], &
transNN(0, 0),
838 int nb_row = data.getIndices().size();
847 invK.resize(3, 3,
false);
852 &
invK(0, 0), &
invK(0, 1), &
invK(0, 2), &
invK(1, 0), &
invK(1, 1),
855 auto t_n_hdiv = data.getFTensor1N<3>();
857 int nb_gauss_pts = data.getN().size1();
858 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
867 for (
int ll = 0; ll != nb_row; ll++) {
868 Nf[ll] +=
w * t_n_hdiv(
i) * t_inv_k(
i,
j) * t_flux(
j);
889 string val_name_col,
Vec f)
904 if (data.getFieldData().size() == 0)
906 int nb_row = data.getIndices().size();
909 divVec.resize(data.getN().size2() / 3, 0);
910 if (
divVec.size() != data.getIndices().size()) {
912 "data inconsistency");
914 int nb_gauss_pts = data.getN().size1();
916 for (; gg < nb_gauss_pts; gg++) {
946 string flux_name_col, Mat aij,
Vec f)
948 val_name_row, flux_name_col,
978 if (
Aij == PETSC_NULL)
980 if (row_data.getFieldData().size() == 0)
982 if (col_data.getFieldData().size() == 0)
984 int nb_row = row_data.getFieldData().size();
985 int nb_col = col_data.getFieldData().size();
986 NN.resize(nb_row, nb_col);
988 divVec.resize(nb_col,
false);
989 int nb_gauss_pts = row_data.getN().size1();
990 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
998 noalias(
NN) +=
w * outer_prod(row_data.getN(gg),
divVec);
1001 &col_data.getIndices()[0], &
NN(0, 0), ADD_VALUES);
1002 transNN.resize(nb_col, nb_row);
1005 &row_data.getIndices()[0], &
transNN(0, 0),
1013 if (data.getIndices().size() == 0)
1015 if (data.getIndices().size() != data.getN().size2()) {
1017 "data inconsistency");
1019 int nb_row = data.getIndices().size();
1022 int nb_gauss_pts = data.getN().size1();
1023 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1053 if (data.getFieldData().size() == 0)
1056 int nb_row = data.getFieldData().size();
1057 Nf.resize(nb_row,
false);
1059 int nb_gauss_pts = data.getN().size1();
1060 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1078 noalias(
Nf) +=
w * data.getN(gg) * flux;
1121 if (data.getFieldData().size() == 0)
1124 nF.resize(data.getIndices().size());
1126 int nb_gauss_pts = data.getN().size1();
1127 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1147 noalias(
nF) +=
w * prod(data.getVectorN<3>(gg),
getNormal()) * value;
1152 &
nF[0], ADD_VALUES);
1186 if (data.getFieldData().size() == 0)
1188 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1189 int nb_dofs = data.getFieldData().size();
1190 int nb_gauss_pts = data.getN().size1();
1191 if (3 * nb_dofs !=
static_cast<int>(data.getN().size2())) {
1193 "wrong number of dofs");
1195 NN.resize(nb_dofs, nb_dofs);
1204 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1206 normal_ptr = &getNormalsAtGaussPts(0)[0];
1209 normal_ptr = &getNormal()[0];
1215 auto t_n_hdiv_row = data.getFTensor1N<3>();
1220 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1224 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1225 x = getHoCoordsAtGaussPts()(gg, 0);
1226 y = getHoCoordsAtGaussPts()(gg, 1);
1227 z = getHoCoordsAtGaussPts()(gg, 2);
1229 x = getCoordsAtGaussPts()(gg, 0);
1230 y = getCoordsAtGaussPts()(gg, 1);
1231 z = getCoordsAtGaussPts()(gg, 2);
1239 double w = getGaussPts()(2, gg);
1241 nrm2 = sqrt(t_normal(
i) * t_normal(
i));
1246 for (
int ll = 0; ll != nb_dofs; ll++) {
1249 &data.getVectorN<3>(gg)(0,
HVEC0),
1250 &data.getVectorN<3>(gg)(0,
HVEC1),
1251 &data.getVectorN<3>(gg)(0,
HVEC2), 3);
1252 for (
int kk = 0; kk <= ll; kk++) {
1253 NN(ll, kk) +=
w * t_n_hdiv_row(
i) * t_n_hdiv_col(
i);
1257 Nf[ll] +=
w * t_n_hdiv_row(
i) * t_normal(
i) * flux / nrm2;
1262 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1264 nrm2 = sqrt(t_normal(
i) * t_normal(
i));
1268 cTx.
bcIndices.insert(data.getIndices().begin(), data.getIndices().end());
1275 CHKERR VecSetOption(
cTx.
D0, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1277 &*
Nf.begin(), INSERT_VALUES);
1278 for (
int dd = 0;
dd != nb_dofs; ++
dd)
1279 data.getFieldDofs()[
dd]->getFieldData() =
Nf[
dd];
1294 const std::string val_name =
"VALUES")
1304 if (data.getFieldData().size() == 0)
1307 int nb_gauss_pts = data.getN().size1();
1309 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1311 inner_prod(trans(data.getN(gg)), data.getFieldData());
1327 const std::string val_name =
"VALUES")
1336 if (data.getFieldData().size() == 0)
1338 int nb_gauss_pts = data.getDiffN().size1();
1340 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1341 ublas::matrix_column<MatrixDouble> values_grad_at_gauss_pts(
1343 noalias(values_grad_at_gauss_pts) =
1344 prod(trans(data.getDiffN(gg)), data.getFieldData());
1359 const std::string field_name)
1369 if (data.getFieldData().size() == 0)
1371 int nb_gauss_pts = data.getDiffN().size1();
1372 int nb_dofs = data.getFieldData().size();
1375 if (
type == MBTRI && side == 0) {
1380 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1383 ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1386 prod(trans(data.getVectorN<3>(gg)), data.getFieldData());
1418 invK.resize(3, 3,
false);
1419 int nb_gauss_pts = data.getN().size1();
1422 Tag th_error_flux, th_error_div;
1424 "ERROR_FLUX", 1, MB_TYPE_DOUBLE, th_error_flux,
1425 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1426 double *error_flux_ptr;
1428 th_error_flux, &fe_ent, 1, (
const void **)&error_flux_ptr);
1431 "ERROR_DIV", 1, MB_TYPE_DOUBLE, th_error_div,
1432 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1433 double *error_div_ptr;
1435 th_error_div, &fe_ent, 1, (
const void **)&error_div_ptr);
1439 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1440 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1441 double *error_jump_ptr;
1443 th_error_jump, &fe_ent, 1, (
const void **)&error_jump_ptr);
1444 *error_jump_ptr = 0;
1447 const double h = pow(
getVolume() * 12 / sqrt(2), (
double)1 / 3);
1449 for (
int ff = 0; ff != 4; ff++) {
1452 double *error_face_jump_ptr;
1454 th_error_jump, &face, 1, (
const void **)&error_face_jump_ptr);
1455 *error_face_jump_ptr = (1 / sqrt(h)) * sqrt(*error_face_jump_ptr);
1456 *error_face_jump_ptr = pow(*error_face_jump_ptr, 2);
1457 *error_jump_ptr += *error_face_jump_ptr;
1460 *error_flux_ptr = 0;
1463 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1484 ublas::matrix_column<MatrixDouble> flux_at_gauss_pt(
1486 ublas::matrix_column<MatrixDouble> values_grad_at_gauss_pts(
1489 prod(
invK, flux_at_gauss_pt) + values_grad_at_gauss_pts;
1492 *error_div_ptr = h * sqrt(*error_div_ptr);
1493 *error_div_ptr = pow(*error_div_ptr, 2);
1494 cTx.
errorMap[sqrt(*error_flux_ptr + *error_div_ptr + *error_jump_ptr)] =
1535 if (data.getFieldData().size() == 0)
1537 int nb_gauss_pts = data.getN().size1();
1538 valMap[getFaceSense()].resize(nb_gauss_pts);
1539 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
1540 valMap[getFaceSense()][gg] =
1541 inner_prod(trans(data.getN(gg)), data.getFieldData());
1560 if (
type == MBTRI) {
1561 EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1566 "ERROR_JUMP", 1, MB_TYPE_DOUBLE, th_error_jump,
1567 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1568 double *error_jump_ptr;
1570 th_error_jump, &fe_ent, 1, (
const void **)&error_jump_ptr);
1571 *error_jump_ptr = 0;
1587 int nb_gauss_pts = data.getN().size1();
1590 if (
valMap.size() == 1) {
1591 if (
static_cast<int>(
valMap.begin()->second.size()) != nb_gauss_pts) {
1593 "wrong number of integration points");
1595 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
1597 if (
static_cast<int>(getNormalsAtGaussPts().size1()) ==
1599 x = getHoCoordsAtGaussPts()(gg, 0);
1600 y = getHoCoordsAtGaussPts()(gg, 1);
1601 z = getHoCoordsAtGaussPts()(gg, 2);
1603 x = getCoordsAtGaussPts()(gg, 0);
1604 y = getCoordsAtGaussPts()(gg, 1);
1605 z = getCoordsAtGaussPts()(gg, 2);
1610 double w = getGaussPts()(2, gg);
1611 if (
static_cast<int>(getNormalsAtGaussPts().size1()) ==
1613 w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1617 *error_jump_ptr +=
w * pow(value -
valMap.begin()->second[gg], 2);
1619 }
else if (
valMap.size() == 2) {
1620 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
1621 double w = getGaussPts()(2, gg);
1622 if (getNormalsAtGaussPts().size1() == (
unsigned int)nb_gauss_pts) {
1623 w *= norm_2(getNormalsAtGaussPts(gg)) * 0.5;
1628 *error_jump_ptr +=
w * pow(
delta, 2);
1632 "data inconsistency, wrong number of neighbors "
1633 "valMap.size() = %d",
1645 #endif //_MIX_TRANPORT_ELEMENT_HPP_