Calculate volume by integrating volume elements and using divergence theorem by integrating surface elements.
static char help[] =
"...\n\n";
Vec vOl;
OpVolume(
const std::string &field_name, Vec vol)
vOl(vol) {}
if (type != MBVERTEX)
const int nb_int_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_ho_det = getFTenosr0HoMeasure();
double v = getMeasure();
double vol = 0;
for (int gg = 0; gg != nb_int_pts; gg++) {
vol += t_w * t_ho_det * v;
++t_w;
++t_ho_det;
}
CHKERR VecSetValue(vOl, 0, vol, ADD_VALUES);
}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
}
};
Vec vOl;
OpFace(
const std::string &field_name, Vec vol)
vOl(vol) {}
if (type != MBVERTEX)
const int nb_int_pts = getGaussPts().size2();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1HoCoordsAtGaussPts();
double vol = 0;
for (int gg = 0; gg != nb_int_pts; gg++) {
vol += (t_coords(i) * t_normal(i)) * t_w;
++t_normal;
++t_w;
++t_coords;
}
vol /= 6;
CHKERR VecSetValue(vOl, 0, vol, ADD_VALUES);
}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
}
};
int operator()(int, int, int) const { return 2; }
};
int operator()(int, int, int) const { return 4; }
};
int main(
int argc,
char *argv[]) {
try {
moab::Core moab_core;
moab::Interface &
moab = moab_core;
DMType dm_name = "DMMOFEM";
{
DM dm;
boost::shared_ptr<ForcesAndSourcesCore> domain_fe =
boost::shared_ptr<ForcesAndSourcesCore>(
boost::shared_ptr<ForcesAndSourcesCore> boundary_fe =
boost::shared_ptr<ForcesAndSourcesCore>(
domain_fe->getRuleHook =
VolRule();
int ghosts[] = {0};
Vec vol, surf_vol;
CHKERR VecCreateGhost(PETSC_COMM_WORLD,
CHKERR VecDuplicate(vol, &surf_vol);
domain_fe->getOpPtrVector().push_back(
new OpVolume(
"MESH_NODE_POSITIONS", vol));
boundary_fe->getOpPtrVector().push_back(
new OpFace(
"MESH_NODE_POSITIONS", surf_vol));
domain_fe);
boundary_fe);
CHKERR VecGhostUpdateBegin(vol, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(vol, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecAssemblyBegin(surf_vol);
CHKERR VecAssemblyEnd(surf_vol);
CHKERR VecGhostUpdateBegin(surf_vol, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(surf_vol, ADD_VALUES, SCATTER_REVERSE);
double *a_vol;
CHKERR VecGetArray(vol, &a_vol);
double *a_surf_vol;
CHKERR VecGetArray(surf_vol, &a_surf_vol);
cout << "Volume = " << a_vol[0] << endl;
cout << "Surf Volume = " << a_surf_vol[0] << endl;
if (fabs(a_vol[0] - a_surf_vol[0]) > 1e-12) {
}
CHKERR VecRestoreArray(vol, &a_vol);
CHKERR VecRestoreArray(vol, &a_surf_vol);
}
}
}
return 0;
}