38 chrono::high_resolution_clock::time_point
start;
39 chrono::high_resolution_clock::time_point
stop;
42 stop = chrono::high_resolution_clock::now();
43 auto duration = chrono::duration_cast<chrono::microseconds>(
stop -
start);
44 cout <<
"Time taken by function: " << duration.count() <<
" ms." << endl;
51 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
53 commonDataPtr(common_data_ptr) {}
58 const size_t nb_dofs = data.
getIndices().size();
62 const size_t nb_base_functions = data.
getN().size2();
63 if (3 * nb_base_functions < nb_dofs)
66 "Number of DOFs is larger than number of base functions on entity");
68 const size_t nb_gauss_pts = data.
getN().size1();
70 auto t_w = getFTensor0IntegrationWeight();
72 auto t_grad = getFTensor2FromMat<3, 3>(*(
commonDataPtr->mGradPtr));
73 auto t_coords = getFTensor1CoordsAtGaussPts();
75 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
78 omg(
i) = (*cache).Omega(
i,
j) * t_coords(
j);
80 double alpha = getMeasure() * t_w * (*cache).density;
81 auto t_nf = getFTensor1FromArray<3, 3>(locF);
84 for (; bb != nb_dofs / 3; ++bb) {
87 t_diff_base(
k) * omg(
k);
92 for (; bb < nb_base_functions; ++bb)
105 const std::string row_field_name,
const std::string col_field_name,
106 boost::shared_ptr<CommonData> common_data_ptr)
109 commonDataPtr(common_data_ptr) {
118 const size_t nb_row_dofs = row_data.
getIndices().size();
119 const size_t nb_col_dofs = col_data.
getIndices().size();
121 if (nb_row_dofs && nb_col_dofs) {
123 const size_t nb_integration_pts = row_data.
getN().size1();
124 const size_t nb_row_base_funcs = row_data.
getN().size2();
126 auto t_w = getFTensor0IntegrationWeight();
127 auto t_coords = getFTensor1CoordsAtGaussPts();
128 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
129 double alpha = getMeasure() * t_w * (*cache).density;
132 omg(
i) = (*cache).Omega(
i,
j) * t_coords(
j);
135 for (; rr != nb_row_dofs / 3; ++rr) {
137 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
141 for (
size_t cc = 0; cc != nb_col_dofs / 3; ++cc) {
142 t_mat(
i,
k) -= alpha * t_row_diff_base(
j) * omg(
j) *
151 for (; rr != nb_row_base_funcs; ++rr)
163 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
164 boost::shared_ptr<DomainSideEle> side_op_fe)
166 commonDataPtr(common_data_ptr), sideOpFe(side_op_fe) {}
171 const size_t nb_gauss_pts = getGaussPts().size2();
172 const size_t nb_dofs = data.
getIndices().size();
176 auto t_w = getFTensor0IntegrationWeight();
177 auto t_coords = getFTensor1CoordsAtGaussPts();
179 auto t_grad = getFTensor2FromMat<3, 3>(*(
commonDataPtr->mGradPtr));
180 size_t nb_base_functions = data.
getN().size2();
182 auto t_normal = getFTensor1Normal();
183 const double ls = sqrt(t_normal(
i) * t_normal(
i));
186 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
188 auto t_nf = getFTensor1FromArray<3, 3>(locF);
190 if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
192 auto t_n = getFTensor1FromPtr<3>(&*
n.data().begin());
193 t_normal(
i) = t_n(
i) / sqrt(t_n(
j) * t_n(
j));
196 const double alpha = t_w * getMeasure() * (*cache).density;
198 omg(
i) = (*cache).Omega(
i,
j) * t_coords(
j);
201 for (; bb != nb_dofs / 3; ++bb) {
203 t_nf(
i) += alpha * t_base *
205 (omg(
k) * t_normal(
k));
211 for (; bb < nb_base_functions; ++bb)
224 const std::string row_field_name,
const std::string col_field_name,
225 boost::shared_ptr<CommonData> common_data_ptr,
226 boost::shared_ptr<DomainSideEle> side_op_fe)
228 commonDataPtr(common_data_ptr), sideOpFe(side_op_fe) {
240 const size_t row_nb_dofs = row_data.
getIndices().size();
241 size_t col_nb_dofs = col_data.
getIndices().size();
243 if (row_nb_dofs && col_nb_dofs) {
245 if (col_type == MBVERTEX) {
256 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
258 commonDataPtr(common_data_ptr) {
272 if (row_type != MBVERTEX)
275 const size_t nb_gauss_pts = getGaussPts().size2();
276 const size_t row_nb_dofs =
commonDataPtr->faceRowData->getIndices().size();
277 const size_t col_nb_dofs = col_data.
getIndices().size();
279 if (row_nb_dofs && col_nb_dofs) {
281 auto t_coords = getFTensor1CoordsAtGaussPts();
282 auto t_w = getFTensor0IntegrationWeight();
283 auto t_row_base =
commonDataPtr->faceRowData->getFTensor0N();
285 auto t_row_diff_base =
commonDataPtr->faceRowData->getFTensor1DiffN<3>();
286 size_t nb_face_functions =
commonDataPtr->faceRowData->getN().size2();
287 auto t_normal = getFTensor1Normal();
288 const double ls = sqrt(t_normal(
i) * t_normal(
i));
291 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
294 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
297 const double alpha = t_w *
commonDataPtr->mMeasure * (*cache).density;
300 omg(
i) = (*cache).Omega(
i,
j) * t_coords(
j);
302 if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
304 auto t_n = getFTensor1FromPtr<3>(&*
n.data().begin());
305 t_normal(
i) = t_n(
i) / sqrt(t_n(
j) * t_n(
j));
309 for (; rr != row_nb_dofs / 3; ++rr) {
310 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
312 const double beta = alpha * t_row_base;
315 for (
size_t cc = 0; cc != col_nb_dofs / 3; ++cc) {
318 (omg(
j) * t_col_diff_base(
j)) * (omg(
l) * t_normal(
l));
326 for (; rr < nb_face_functions; ++rr)
339 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
341 commonDataPtr(common_data_ptr) {}
347 const size_t nb_dofs = data.
getIndices().size();
349 const size_t nb_integration_pts = data.
getN().size1();
351 auto t_plastic_strain_dot =
352 getFTensor2SymmetricFromMat<3>(*(
commonDataPtr->plasticStrainDotPtr));
359 auto t_grad_plastic_strain =
360 getFTensor3DgFromMat<3, 3>(*(
commonDataPtr->plasticGradStrainPtr));
361 commonDataPtr->guidingVelocityPtr->resize(3, nb_integration_pts,
false);
363 auto t_omega = getFTensor1FromMat<3>(*
commonDataPtr->guidingVelocityPtr);
365 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
367 t_omega(
i) = (*cache).Omega(
i,
j) * t_coords(
j);
368 t_tau_dot += t_grad_tau(
i) * t_omega(
i);
369 t_plastic_strain_dot(
i,
j) += t_grad_plastic_strain(
i,
j,
k) * t_omega(
k);
373 ++t_grad_plastic_strain;
374 ++t_plastic_strain_dot;
384 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
385 boost::shared_ptr<DomainSideEle> side_op_fe)
387 commonDataPtr(common_data_ptr), sideOpFe(side_op_fe) {
388 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
389 doEntities[MBTRI] = doEntities[MBQUAD] =
true;
397 if ((type == MBTRI || type == MBQUAD) && side == 0) {
398 const size_t nb_dofs = data.
getIndices().size();
402 const size_t nb_integration_pts = data.
getN().size1();
403 const size_t nb_base_functions = data.
getN().size2();
406 const int check_size =
commonDataPtr->plasticStrainSideMap.size();
409 if (
static_cast<int>(
410 commonDataPtr->plasticStrainSideMap.begin()->second.size2()) !=
413 "wrong number of integration points");
415 auto &plastic_mat_l =
commonDataPtr->plasticStrainSideMap.at(-1);
416 auto &plastic_mat_r =
commonDataPtr->plasticStrainSideMap.at(1);
417 auto t_plastic_strain_l = getFTensor2SymmetricFromMat<3>(plastic_mat_l);
418 auto t_plastic_strain_r = getFTensor2SymmetricFromMat<3>(plastic_mat_r);
424 auto t_w = getFTensor0IntegrationWeight();
426 auto t_coords = getFTensor1CoordsAtGaussPts();
427 commonDataPtr->plasticTauJumpPtr->resize(nb_integration_pts,
false);
428 commonDataPtr->plasticStrainJumpPtr->resize(6, nb_integration_pts,
false);
429 auto t_plastic_strain_jump =
430 getFTensor2SymmetricFromMat<3>(*
commonDataPtr->plasticStrainJumpPtr);
433 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
434 double alpha = getMeasure() * t_w;
441 t_plastic_strain_jump(
i,
j) =
442 t_plastic_strain_l(
i,
j) - t_plastic_strain_r(
i,
j);
443 t_tau_jump = t_tau_l - t_tau_r;
455 ++t_plastic_strain_jump;
457 ++t_plastic_strain_l;
458 ++t_plastic_strain_r;
471 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
473 commonDataPtr(common_data_ptr) {
474 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
475 doEntities[MBTET] =
true;
476 doEntities[MBHEX] =
true;
482 const int nb_gauss_pts = getGaussPts().size2();
483 const int nb_base_functions = data.
getN().size2();
490 const size_t nb_dofs = data.
getIndices().size();
496 MatrixDouble diff = getCoordsAtGaussPts() - getFaceCoordsAtGaussPts();
497 const double eps = 1e-12;
498 if (norm_inf(diff) >
eps)
500 "coordinates at integration pts are different");
502 const size_t nb_gauss_pts = data.
getN().size1();
506 int test_nb_in_loop = getFEMethod()->nInTheLoop;
507 int test_face_sense = getFaceSense();
509 mat_ep.resize(6, nb_gauss_pts,
false);
519 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<3>(mat_ep);
520 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
523 for (; bb != nb_dofs / 6; ++bb) {
524 values_at_gauss_pts(
i,
j) += field_data(
i,
j) * base_function;
528 for (; bb != nb_base_functions; ++bb)
530 ++values_at_gauss_pts;
537 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr)
539 commonDataPtr(common_data_ptr) {
540 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
541 doEntities[MBTET] =
true;
542 doEntities[MBHEX] =
true;
550 const size_t nb_dofs = data.
getIndices().size();
553 int nb_gauss_pts = data.
getN().size1();
555 vec_tau.resize(nb_gauss_pts,
false);
558 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
566 const std::string
field_name, moab::Interface &post_proc_mesh,
567 std::vector<EntityHandle> &map_gauss_pts,
568 boost::shared_ptr<CommonData> common_data_ptr)
570 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
571 commonDataPtr(common_data_ptr) {
573 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
574 doEntities[MBTRI] = doEntities[MBQUAD] =
true;
583 const size_t nb_gauss_pts = getGaussPts().size2();
584 std::array<double, 9> def;
585 std::fill(def.begin(), def.end(), 0);
587 auto get_tag = [&](
const std::string name,
size_t size) {
590 MB_TAG_CREAT | MB_TAG_SPARSE,
599 for (
size_t r = 0; r != 3; ++r)
600 for (
size_t c = 0;
c != 3; ++
c)
606 for (
size_t r = 0; r != 3; ++r)
618 &*mat.data().begin());
621 auto th_rot_vel = get_tag(
"ROT_VELOCITY", 3);
622 auto th_tau_jump = get_tag(
"TAU_JUMP", 1);
623 auto th_plastic_jump = get_tag(
"PLASTIC_STRAIN_JUMP", 9);
624 auto t_coords = getFTensor1CoordsAtGaussPts();
626 getFTensor2SymmetricFromMat<3>(*(
commonDataPtr->plasticStrainJumpPtr));
629 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
632 velocity(
i) = (*cache).Omega(
i,
j) * t_coords(
j);
634 CHKERR set_tag(th_tau_jump, gg, set_scalar(t_tau_jump));
635 CHKERR set_tag(th_plastic_jump, gg, set_matrix_3d(t_ep_jump));
636 CHKERR set_tag(th_rot_vel, gg, set_vector(velocity));
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
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
Tensors class implemented by Walter Landry.
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
constexpr double t
plate stiffness
constexpr auto field_name
chrono::high_resolution_clock::time_point stop
chrono::high_resolution_clock::time_point start
bool sYmm
If true assume that matrix is symmetric structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(Tensor_Dim *(Tensor_Dim+1))/2 >, Tensor_Dim > getFTensor2SymmetricFieldData()
Return symmetric FTensor rank 2, i.e. matrix from filed data coeffinects.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
default operator for TRI element
MoFEMErrorCode loopSideVolumes(const string fe_name, VolumeElementForcesAndSourcesCoreOnSide &fe_method)
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
double getMeasure() const
get measure of element
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
default operator for TET element
OpPostProcPlasticSkeleton(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
boost::shared_ptr< CommonData > commonDataPtr
std::vector< EntityHandle > & mapGaussPts
moab::Interface & postProcMesh