78 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Topo Config",
"none");
84 CHKERR PetscOptionsScalar(
"-young",
"default yOung modulus",
"",
yOung, &
yOung,
87 CHKERR PetscOptionsScalar(
"-poisson",
"poisson ratio",
"", 0.1, &
pOisson,
90 CHKERR PetscOptionsScalar(
"-move",
"max move",
"", 0.1, &
mOve, PETSC_NULL);
93 double lambda_coef[3] = {0.,0.,0.};
95 PetscBool is_lambda_defined = PETSC_FALSE;
96 PetscOptionsGetRealArray(NULL,
"",
"-lambda", lambda_coef,
97 &nmax, &is_lambda_defined);
114 "Wrong number of lambda parameters, should be 1 or 3, data "
123 ierr = PetscOptionsEnd();
130 D.resize(36, 1,
false);
140 t_D(
i,
j,
k,
l) = 0.;
146 t_D(0, 1, 0, 1) = 0.5 * (1 - 2 *
pOisson);
147 t_D(0, 2, 0, 2) = 0.5 * (1 - 2 *
pOisson);
148 t_D(1, 2, 1, 2) = 0.5 * (1 - 2 *
pOisson);
157 t_D(
i,
j,
k,
l) *= coefficient;
168 thSensitivity, &fe_ent, 1, (
const void **)&tag_data, &tag_size);
176 tag_data, &tag_size);
177 }
else if (tag_size != nb_gauss_pts) {
179 "Wrong size of the tag, data inconsistency");
182 tag_size, ublas::shallow_array_adaptor<double>(tag_size, tag_data));
194 tag_data, &tag_size);
203 thDens, &fe_ent, 1, (
const void **)&tag_data, &tag_size);
208 void const *tag_data[] = {&*
rhoVecAtTag.data().begin()};
211 tag_data, &tag_size);
212 }
else if (tag_size != nb_gauss_pts) {
214 "Wrong size of the tag, data inconsistency");
217 tag_size, ublas::shallow_array_adaptor<double>(tag_size, tag_data));
225 void const *tag_data[] = {&*
rhoVecAtTag.data().begin()};
228 tag_data, &tag_size);
234 std::vector<double> def_val(9, 0);
237 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, &def_val[0]);
239 "RHO", 1, MB_TYPE_DOUBLE,
thDens,
240 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, &def_val[0]);
246 PetscBool is_partitioned) {
248 DMType dm_name =
"DM_TOPO";
250 CHKERR DMRegister_MoFEM(dm_name);
251 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
252 CHKERR DMSetType(dm, dm_name);
253 CHKERR DMMoFEMCreateMoFEM(dm, &m_field, dm_name, bit_level0);
254 CHKERR DMSetFromOptions(dm);
255 CHKERR DMMoFEMAddElement(dm,
"ELASTIC");
256 CHKERR DMMoFEMAddElement(dm,
"RHO_APPROX");
258 CHKERR DMMoFEMAddElement(dm,
"PRESSURE_FE");
260 CHKERR DMMoFEMAddElement(dm,
"FORCE_FE");
261 CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
273 PetscBool is_partitioned) {
276 DMType dm_name_sub_elast =
"DM_ELASTIC";
277 CHKERR DMRegister_MoFEM(dm_name_sub_elast);
278 CHKERR DMCreate(PETSC_COMM_WORLD, &sub_dm_elastic);
279 CHKERR DMSetType(sub_dm_elastic, dm_name_sub_elast);
282 CHKERR DMMoFEMCreateSubDM(sub_dm_elastic, dm,
"ELASTIC_PROB");
283 CHKERR DMMoFEMSetDestroyProblem(sub_dm_elastic, PETSC_TRUE);
284 CHKERR DMMoFEMAddSubFieldRow(sub_dm_elastic,
"DISPLACEMENTS");
285 CHKERR DMMoFEMAddSubFieldCol(sub_dm_elastic,
"DISPLACEMENTS");
286 CHKERR DMMoFEMAddElement(sub_dm_elastic,
"ELASTIC");
288 CHKERR DMMoFEMAddElement(sub_dm_elastic,
"PRESSURE_FE");
290 CHKERR DMMoFEMAddElement(sub_dm_elastic,
"FORCE_FE");
292 CHKERR DMMoFEMSetSquareProblem(sub_dm_elastic, PETSC_TRUE);
293 CHKERR DMMoFEMSetIsPartitioned(sub_dm_elastic, is_partitioned);
294 CHKERR DMSetUp(sub_dm_elastic);
296 DMType dm_name_sub_rho =
"DM_DENSITY";
297 CHKERR DMRegister_MoFEM(dm_name_sub_rho);
298 CHKERR DMCreate(PETSC_COMM_WORLD, &sub_dm_rho);
299 CHKERR DMSetType(sub_dm_rho, dm_name_sub_rho);
302 CHKERR DMMoFEMCreateSubDM(sub_dm_rho, dm,
"DENSITY_PROB");
303 CHKERR DMMoFEMSetDestroyProblem(sub_dm_rho, PETSC_TRUE);
304 CHKERR DMMoFEMAddSubFieldRow(sub_dm_rho,
"RHO");
305 CHKERR DMMoFEMAddSubFieldCol(sub_dm_rho,
"RHO");
306 CHKERR DMMoFEMAddElement(sub_dm_rho,
"RHO_APPROX");
307 CHKERR DMMoFEMSetSquareProblem(sub_dm_rho, PETSC_TRUE);
308 CHKERR DMMoFEMSetIsPartitioned(sub_dm_rho, is_partitioned);
309 CHKERR DMSetUp(sub_dm_rho);
314 std::vector<double> &coords,
double &boxSize) {
317 material_coords[0], material_coords[1], material_coords[2]);
321 double max_xyz[] = {0, 0, 0};
322 for (
int bb = 0; bb != coords.size() / 3; bb++) {
323 t_coords(
i) -= t_material_coords(
i);
324 for (
int dd = 0; dd != 3; dd++) {
325 max_xyz[dd] = (max_xyz[dd] < fabs(t_coords(dd))) ? fabs(t_coords(dd))
331 boxSize = (max_xyz[0] > max_xyz[1]) ? max_xyz[0] : max_xyz[1];
332 boxSize = (boxSize > max_xyz[2]) ? boxSize : max_xyz[2];
339 :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
350 "DISPLACEMENTS",
"DISPLACEMENTS", OPROWCOL, symm),
365 DataForcesAndSourcesCore::EntData &row_data,
366 DataForcesAndSourcesCore::EntData &col_data) {
369 nbRows = row_data.getIndices().size();
374 nbCols = col_data.getIndices().size();
381 if (row_side == col_side && row_type == col_type) {
409 virtual MoFEMErrorCode
411 DataForcesAndSourcesCore::EntData &col_data) {
414 int nb_dofs_row = row_data.getFieldData().size();
415 if (nb_dofs_row == 0)
417 int nb_dofs_col = col_data.getFieldData().size();
418 if (nb_dofs_col == 0)
421 K.resize(nb_dofs_row, nb_dofs_col,
false);
423 auto max = [](
double a,
double b) {
return a > b ?
a : b; };
424 auto min = [](
double a,
double b) {
return a < b ?
a : b; };
427 auto get_tensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
429 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
430 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
431 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2));
440 double vol = getVolume();
443 auto t_w = getFTensor0IntegrationWeight();
445 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
446 const int row_nb_base_fun = row_data.getN().size2();
454 double a = t_w * vol;
455 if (getHoGaussPtsDetJac().size()) {
456 a *= getHoGaussPtsDetJac()[gg];
460 for (; rr !=
nbRows / 3; ++rr) {
461 auto t_m = get_tensor2(
K, 3 * rr, 0);
462 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
464 t_rowD(
l,
j,
k) = t_D(
i,
j,
k,
l) * (
a * t_row_diff_base(
i));
465 for (
int cc = 0; cc !=
nbCols / 3; ++cc) {
466 t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base(
k);
472 for (; rr != row_nb_base_fun; ++rr)
486 virtual MoFEMErrorCode
aSsemble(DataForcesAndSourcesCore::EntData &row_data,
487 DataForcesAndSourcesCore::EntData &col_data) {
490 int nb_dofs_row = row_data.getFieldData().size();
491 if (nb_dofs_row == 0)
493 int nb_dofs_col = col_data.getFieldData().size();
494 if (nb_dofs_col == 0)
498 const int *row_indices = &*row_data.getIndices().data().begin();
500 const int *col_indices = &*col_data.getIndices().data().begin();
501 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
502 : getFEMethod()->snes_B;
505 &*
K.data().begin(), ADD_VALUES);
509 &*
K.data().begin(), ADD_VALUES);
527 DataForcesAndSourcesCore::EntData &data) {
529 if (type != MBVERTEX)
531 int nb_dofs = data.getFieldData().size();
556 double newSensitivity = 0;
557 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
560 strain(
i,
j) = 0.5 * (grad(
i,
j) || grad(
j,
i));
561 const double trace = strain(
i,
i);
566 t_sensitivity = (-
p) * pow(
rho,
p - 1) * energy * vol;
570 if (t_sensitivity > 0)
573 "t_sensitivity is not within boundaries, that should not happen");
590template <
bool UPDATE = true>
600 DataForcesAndSourcesCore::EntData &data) {
602 if (type != MBVERTEX)
604 int nb_dofs = data.getFieldData().size();
609 int nb_gauss_pts = data.getN().size1();
620 auto max = [](
double a,
double b) {
return a > b ?
a : b; };
621 auto min = [](
double a,
double b) {
return a < b ?
a : b; };
625 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
627 VectorAdaptor shape_function = data.getN(gg, nb_dofs);
630 const double sensitivity = t_sens;
632 double rho = rho_old ;
634 const double &x =
rho;
638 const double xnew =
rho * sqrt(-sensitivity/den);
640 rho_new = max(0.001, min(1.0, max(x - mv, min(1.0, min(x + mv, xnew)))));
641 if (rho_new < 0.001 || rho_new > 1.00001)
643 "rho_new is not within boundaries, that should not happen");
644 if (rho_old < 0.001 || rho_old > 1.00001)
646 "rho_new is not within boundaries, that should not happen");
647 const double &check = rho_new;
648 if (!std::isnormal(check))
650 "rho_new is not normal, that should not happen");
652 change += vol * fabs(rho_old - rho_new);
676 DataForcesAndSourcesCore::EntData &row_data) {
678 if (row_type != MBVERTEX)
700 DataForcesAndSourcesCore::EntData &row_data) {
702 if (row_type != MBVERTEX)
704 int nb_gauss_pts = row_data.getN().size1();
707 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
710 if (getHoGaussPtsDetJac().size() > 0) {
711 vol *= getHoGaussPtsDetJac()[gg];
739 DataForcesAndSourcesCore::EntData &row_data,
740 DataForcesAndSourcesCore::EntData &col_data) {
744 int nb_rows = row_data.getIndices().size();
745 int nb_cols = col_data.getIndices().size();
751 locK.resize(nb_rows, nb_cols,
false);
765 int nb_gauss_pts = row_data.getN().size1();
766 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
768 if (getHoGaussPtsDetJac().size() > 0) {
769 vol *= getHoGaussPtsDetJac()[gg];
771 VectorAdaptor row_base_function = row_data.getN(gg);
772 VectorAdaptor col_base_function = col_data.getN(gg);
773 locK += vol * outer_prod(row_base_function, col_base_function);
775 MatrixAdaptor row_diff_base_function = row_data.getDiffN(gg);
776 MatrixAdaptor col_diff_base_function = col_data.getDiffN(gg);
778 row_diff_base_function = prod(row_diff_base_function, anisotropic_filter);
779 col_diff_base_function = prod(col_diff_base_function, anisotropic_filter);
782 locK += vol * prod(row_diff_base_function, trans(col_diff_base_function));
785 &row_data.getIndices()[0], col_data.getIndices().size(),
786 &col_data.getIndices()[0], &
locK(0, 0), ADD_VALUES);
789 if (row_side != col_side || row_type != col_type) {
793 &col_data.getIndices()[0],
794 row_data.getIndices().size(),
795 &row_data.getIndices()[0], &
transK(0, 0), ADD_VALUES);
814 DataForcesAndSourcesCore::EntData &row_data) {
817 int nb_rows = row_data.getIndices().size();
820 nF.resize(nb_rows,
false);
825 int nb_gauss_pts = row_data.getN().size1();
827 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
830 if (getHoGaussPtsDetJac().size() > 0) {
831 vol *= getHoGaussPtsDetJac()[gg];
833 VectorAdaptor base_function = row_data.getN(gg);
835 nF += vol *
rho * base_function;
837 if (
rho > 1.000001 ||
rho < 0.001)
839 "wrong rho. data inconsistency");
845 &row_data.getIndices()[0], &
nF[0], ADD_VALUES);
858 std::vector<EntityHandle> &map_gauss_pts,
865 DataForcesAndSourcesCore::EntData &data) {
867 if (type != MBVERTEX)
870 bzero(def_VAL, 9 *
sizeof(
double));
881 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
887 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
889 strain(
i,
j) = 0.5 * (grad(
i,
j) || grad(
j,
i));
890 const double trace = grad(
i,
i);
891 const double energy =
908 Mat aij, Vec
x, Vec
f)
917 ParallelComm *pcomm =
921 string name = it->getName();
922 if (name.compare(0, 7,
"DENSITY") == 0) {
925 VectorDouble scaled_values(1);
926 scaled_values[0] = 1;
935 moab::Interface::UNION);
936 ents.insert(_edges.begin(), _edges.end());
941 ents.insert(_nodes.begin(), _nodes.end());
943 auto &dofs_by_uid =
problemPtr->getNumeredRowDofs()->get<Unique_mi_tag>();
944 for (
auto eit = ents.pair_begin(); eit != ents.pair_end(); ++eit) {
946 auto lo_dit = dofs_by_uid.lower_bound(
947 DofEntity::getLoFieldEntityUId(bit_number, eit->first));
948 auto hi_dit = dofs_by_uid.upper_bound(
949 DofEntity::getHiFieldEntityUId(bit_number, eit->second));
950 for (; lo_dit != hi_dit; ++lo_dit) {
952 if (dof->getEntType() == MBVERTEX) {
953 mapZeroRows[dof->getPetscGlobalDofIdx()] = scaled_values[0];
979 std::map<DofIdx, FieldData>::iterator mit =
mapZeroRows.begin();
1013 dIag, PETSC_NULL, PETSC_NULL);
1018 for (std::vector<int>::iterator vit =
dofsIndices.begin();
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static PetscErrorCode ierr
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ 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 ...
double trace(FTensor::Tensor2< T, 2, 2 > &t_stress)
FTensor::Index< 'm', SPACE_DIM > m
#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 c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorShallowArrayAdaptor< double > VectorAdaptor
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
constexpr auto field_name
boost::shared_ptr< VectorDouble > dataRhotGaussPtr
MoFEMErrorCode build_main_dm(DM &dm, DM &sub_dm_elastic, DM &sub_dm_rho, MoFEM::Interface &m_field, BitRefLevel &bit_level0, PetscBool is_partitioned)
MoFEM::Interface & mField
VectorDouble sEnsitivityVec
CommonData(MoFEM::Interface &m_field)
vector< VectorDouble > sTress
MoFEMErrorCode getParameters()
MoFEMErrorCode setSensitivityToTag(const EntityHandle fe_ent, const double sens)
MatrixDouble anisLambdaMat
MoFEMErrorCode build_sub_dms(DM &dm, DM &sub_dm_elastic, DM &sub_dm_rho, MoFEM::Interface &m_field, BitRefLevel &bit_level0, PetscBool is_partitioned)
MoFEMErrorCode getSensitivityFromTag(const EntityHandle fe_ent, const int nb_gauss_pts)
MoFEMErrorCode getBoundBox(const double material_coords[3], std::vector< double > &coords, double &boxSize)
boost::shared_ptr< MatrixDouble > dispGradAtGaussPts
MatrixDouble cOmplianceMat
MoFEMErrorCode setDensityToTag(const EntityHandle fe_ent, const double rho)
MoFEMErrorCode createTags(moab::Interface &moab)
MoFEMErrorCode getDensityFromTag(const EntityHandle fe_ent, const int nb_gauss_pts)
DirichletDensityBc(MoFEM::Interface &m_field, const std::string &field_name, Mat aij, Vec x, Vec f)
MoFEMErrorCode iNitalize()
MoFEMErrorCode postProcess()
function is run at the end of loop
DirichletDensityBc(MoFEM::Interface &m_field, const std::string &field_name)
Set Dirichlet boundary conditions on displacements.
std::map< DofIdx, FieldData > mapZeroRows
const std::string fieldName
field name to set Dirichlet BC
double dIag
diagonal value set on zeroed column and rows
std::vector< int > dofsIndices
boost::ptr_vector< MethodForForceScaling > methodsOp
std::vector< double > dofsValues
MoFEM::Interface & mField
const Problem * problemPtr
raw pointer to problem
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Problem manager is used to build and partition problems.
Mat & snes_B
preconditioner of jacobian matrix
Vec & ts_F
residual vector
Mat & ts_B
Preconditioner for ts_A.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
double getVolume() const
element volume (linear geometry)
OpAssembleDensityRhs(const string &row_field, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
virtual MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate grad-grad operator.
int nbIntegrationPts
number of integration points
OpAssembleK(CommonData &common_data, bool symm=true)
int nbCols
number if dof on column
bool isDiag
true if this block is on diagonal
virtual MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble local entity block matrix.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Do calculations for give operator.
OpAssembleRhoLhs(const std::string &row_field, const std::string &col_field, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpCalculateSensitivities(CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpMassCalculation(const string &field_name, CommonData &common_data, Vec mass_vec)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpPostProcCompliance(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, CommonData &common_data)
std::vector< EntityHandle > & mapGaussPts
moab::Interface & postProcMesh
OpUpdateDensity(CommonData &common_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpVolumeCalculation(Vec &volume_vec)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)