7 #ifndef __SURFACE_SLIDING_CONSTRAINS_HPP__
8 #define __SURFACE_SLIDING_CONSTRAINS_HPP__
21 struct OpGetActiveDofsLambda
25 boost::shared_ptr<VectorDouble> &active_variables_ptr)
32 if (
type == MBVERTEX) {
33 for (
unsigned int dd = 0;
dd != data.getFieldData().size(); ++
dd)
40 template <
int SizeLambda>
41 struct OpGetActiveDofsPositions
46 boost::shared_ptr<VectorDouble> &active_variables_ptr)
53 if (
type == MBVERTEX) {
54 for (
unsigned int dd = 0;
dd != data.getFieldData().size(); ++
dd)
61 template <
int SizeLambda,
int SizePositions>
67 boost::shared_ptr<VectorDouble> &results_ptr)
81 else if (indices.size() == SizeLambda)
83 else if (indices.size() == SizePositions)
87 "Element %s: Data inconsistency nb of indices %d",
90 CHKERR VecSetOption(
getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
96 template <
int SizeLambda,
int SizePositions>
102 const std::string field_name_col,
103 boost::shared_ptr<MatrixDouble> &jacobian_ptr)
115 if (row_type != MBVERTEX)
117 if (col_type != MBVERTEX)
119 VectorInt &row_indices = row_data.getIndices();
120 VectorInt &col_indices = col_data.getIndices();
121 if (row_indices.empty() || col_indices.empty())
125 if (row_indices.size() == SizeLambda)
127 else if (row_indices.size() == SizePositions)
128 shift_row = SizeLambda;
131 "Data inconsistency");
134 if (col_indices.size() == SizeLambda)
136 else if (col_indices.size() == SizePositions)
137 shift_col = SizeLambda;
140 "Data inconsistency");
142 MatrixDouble jac(row_indices.size(), col_indices.size());
143 for (
unsigned int rr = 0; rr != row_indices.size(); ++rr) {
144 for (
unsigned int cc = 0; cc != col_indices.size(); ++cc) {
145 jac(rr, cc) = (*jacobianPtr)(shift_row + rr, shift_col + cc);
146 if (jac(rr, cc) != jac(rr, cc))
148 "Jacobian assemble not a number jac(rr,cc) = %3.4f",
153 col_indices.size(), &col_indices[0], &jac(0, 0),
334 struct DriverElementOrientation {
370 if (
B != PETSC_NULL) {
374 if (
F != PETSC_NULL) {
394 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
395 "Get surface sliding constrains element scaling",
398 CHKERR PetscOptionsScalar(
"-surface_sliding_alpha",
"scaling parameter",
"",
400 ierr = PetscOptionsEnd();
412 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
430 boost::shared_ptr<VectorDouble> &active_variables_ptr,
431 boost::shared_ptr<VectorDouble> &results_ptr,
432 boost::shared_ptr<MatrixDouble> &jacobian_ptr,
445 if (
type != MBVERTEX)
452 ublas::vector<adouble> lambda_dofs(3);
453 for (
int dd = 0;
dd != 3; ++
dd) {
454 lambda_dofs[
dd] <<= (*activeVariablesPtr)[
dd];
456 ublas::vector<adouble> position_dofs(9);
457 for (
int dd = 0;
dd != 9; ++
dd) {
458 position_dofs[
dd] <<= (*activeVariablesPtr)[3 +
dd];
472 int nb_gauss_pts = data.getN().size1();
473 int nb_base_functions = data.getN().size2();
474 auto t_base1 = data.getFTensor0N();
475 auto t_base2 = data.getFTensor0N();
476 auto t_diff_base = data.getFTensor1DiffN<2>();
484 ublas::vector<adouble> c_vec(3);
485 ublas::vector<adouble> f_vec(9);
491 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
494 &position_dofs[0], &position_dofs[1], &position_dofs[2], 3);
498 t_position_ksi(
i) = 0;
499 t_position_eta(
i) = 0;
502 for (
int bb = 0; bb != nb_base_functions; ++bb) {
503 t_position(
i) += t_base1 * t_position_dofs(
i);
504 t_position_ksi(
i) += t_diff_base(N0) * t_position_dofs(
i);
505 t_position_eta(
i) += t_diff_base(N1) * t_position_dofs(
i);
506 lambda += t_base1 * t_lambda_dof;
513 t_delta(
i) = t_position(
i) - t_coord_ref(
i);
521 adouble area = sqrt(t_normal(
i) * t_normal(
i));
523 for (
int bb = 0; bb != nb_base_functions; ++bb) {
524 if (indices[bb] != -1) {
525 val =
w * eo * t_base2;
526 t_c += val * t_normal(
i) * t_delta(
i);
528 t_f(
i) += val * t_normal(
i);
538 for (
int rr = 0; rr != 3; ++rr)
541 for (
int rr = 0; rr != 9; ++rr)
542 f_vec(rr) >>= (*resultsPtr)[3 + rr];
547 double *jac_ptr[3 + 9];
548 for (
int rr = 0; rr != 3 + 9; ++rr)
556 "ADOL-C function evaluation with error");
564 const std::string lagrange_multipliers_field_name,
565 const std::string material_field_name,
566 const double *alpha =
nullptr) {
569 if (alpha !=
nullptr) {
573 boost::shared_ptr<VectorDouble> active_variables_ptr(
575 boost::shared_ptr<VectorDouble> results_ptr(
new VectorDouble(3 + 9));
576 boost::shared_ptr<MatrixDouble> jacobian_ptr(
581 lagrange_multipliers_field_name, active_variables_ptr));
583 material_field_name, active_variables_ptr));
585 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
595 lagrange_multipliers_field_name, active_variables_ptr));
597 material_field_name, active_variables_ptr));
599 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
602 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
604 material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
606 material_field_name, material_field_name, jacobian_ptr));
613 const std::string lagrange_multipliers_field_name,
614 const std::string material_field_name) {
617 boost::shared_ptr<VectorDouble> active_variables_ptr(
619 boost::shared_ptr<VectorDouble> results_ptr(
new VectorDouble(3 + 9));
620 boost::shared_ptr<MatrixDouble> jacobian_ptr(
626 lagrange_multipliers_field_name, active_variables_ptr));
628 material_field_name, active_variables_ptr));
630 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
633 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
641 struct CalculateEdgeBase {
644 Tag &th2, Tag &th3) {
646 double def_val[] = {0, 0, 0};
647 CHKERR moab.tag_get_handle(
"EDGE_BASE0", 3, MB_TYPE_DOUBLE, th0,
648 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
649 CHKERR moab.tag_get_handle(
"EDGE_BASE1", 3, MB_TYPE_DOUBLE, th1,
650 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
651 int def_numb_val[] = {-1};
652 CHKERR moab.tag_get_handle(
"PATCH_NUMBER", 1, MB_TYPE_INTEGER, th2,
653 MB_TAG_CREAT | MB_TAG_SPARSE, def_numb_val);
654 int def_orientation_val[] = {1};
655 CHKERR moab.tag_get_handle(
"PATCH_ORIENTATION", 1, MB_TYPE_INTEGER, th3,
656 MB_TAG_CREAT | MB_TAG_SPARSE,
657 def_orientation_val);
666 auto get_edges = [&](
const Range &ents) {
668 CHKERR moab.get_adjacencies(ents, 1,
false, edges,
669 moab::Interface::UNION);
673 auto get_face_adj = [&, get_edges](
const Range &faces) {
675 CHKERR moab.get_adjacencies(subtract(get_edges(faces), edges), 2,
false,
676 adj_faces, moab::Interface::UNION);
677 return intersect(adj_faces, tris);
680 auto get_patch = [&, get_face_adj](
const EntityHandle face) {
682 patch_ents.insert(face);
685 nb0 = patch_ents.size();
686 patch_ents.merge(get_face_adj(patch_ents));
687 }
while (nb0 != patch_ents.size());
691 auto get_patches = [&]() {
692 std::vector<Range> patches;
693 while (!tris.empty()) {
694 patches.push_back(get_patch(tris[0]));
695 tris = subtract(tris, patches.back());
700 Tag th0, th1, th2, th3;
703 auto patches = get_patches();
705 for (
auto patch : patches) {
708 std::vector<int> tags_vals(patch.size(), pp);
709 CHKERR moab.tag_set_data(th2, patch, &*tags_vals.begin());
718 bool number_pathes =
true,
719 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
720 surface_orientation =
nullptr,
723 Tag th0, th1, th2, th3;
728 for (
auto edge : edges) {
730 CHKERR moab.get_adjacencies(&edge, 1, 2,
false, adj_faces);
731 adj_faces = intersect(adj_faces, tris);
732 if (adj_faces.size() != 2) {
734 "Should be 2 faces adjacent to edge but is %d",
740 &
v[0], &
v[1], &
v[2]);
747 std::array<double, 3> areas;
748 auto calculate_normals = [&, get_tensor_from_vec]() {
750 for (
auto face : adj_faces) {
751 double &x = (
v[ff])[0];
752 double &y = (
v[ff])[1];
753 double &z = (
v[ff])[2];
754 moab::Util::normal(&moab, face, x, y, z);
755 auto t_n = get_tensor_from_vec(
v[ff]);
756 areas[ff] = sqrt(t_n(
i) * t_n(
i));
760 CHKERR moab.tag_get_data(th3, &face, 1, &orientation);
761 if (orientation == -1) {
764 if (surface_orientation && m_field_ptr) {
765 CHKERR surface_orientation->getElementOrientation(*m_field_ptr,
767 int eo = surface_orientation->elementOrientation;
777 auto get_patch_number = [&]() {
778 std::vector<int> p = {0, 0};
779 CHKERR moab.tag_get_data(th2, adj_faces, &*p.begin());
783 std::array<EntityHandle, 2> faces_handles = {adj_faces[0],
785 auto order_normals = [&, get_patch_number]() {
786 auto p = get_patch_number();
789 faces_handles[0] = adj_faces[1];
790 faces_handles[1] = adj_faces[0];
796 auto t_n0 = get_tensor_from_vec(
v[0]);
797 auto t_n1 = get_tensor_from_vec(
v[1]);
800 auto get_base_for_coplanar_faces = [&]() {
803 std::vector<EntityHandle> face_conn;
804 CHKERR moab.get_connectivity(&faces_handles[1], 1, face_conn,
true);
805 std::vector<EntityHandle> edge_conn;
806 CHKERR moab.get_connectivity(&edge, 1, edge_conn,
true);
807 std::sort(face_conn.begin(), face_conn.end());
808 std::sort(edge_conn.begin(), edge_conn.end());
811 if (face_conn[
n] != edge_conn[
n])
814 CHKERR moab.get_coords(&face_conn[
n], 1,
815 &*coords_face.data().begin());
817 CHKERR moab.get_coords(&edge_conn[0], 2,
818 &*coords_edge.data().begin());
820 coords_edge[0], coords_edge[1], coords_edge[2]);
822 coords_edge[3], coords_edge[4], coords_edge[5]);
823 auto t_face_n = get_tensor_from_vec(coords_face);
824 t_face_n(
i) -= t_edge_n0(
i);
826 t_delta(
i) = t_edge_n1(
i) - t_edge_n0(
i);
827 t_delta(
i) /= sqrt(t_delta(
i) * t_delta(
i));
829 if (t_n1(
i) * t_face_n(
i) < 0)
835 auto get_base_for_not_planar_faces = [&]() {
838 t_n1(
i) /= sqrt(t_n1(
i) * t_n1(
i));
843 constexpr
double tol = 1e-6;
844 if ((t_cross(
i) * t_cross(
i)) <
tol)
845 CHKERR get_base_for_coplanar_faces();
847 CHKERR get_base_for_not_planar_faces();
851 CHKERR moab.tag_set_data(th0, &edge, 1, &v0[0]);
852 CHKERR moab.tag_set_data(th1, &edge, 1, &v1[0]);
863 CHKERR moab.create_meshset(MESHSET_SET, meshset);
864 CHKERR moab.add_entities(meshset, edges);
865 if (faces !=
nullptr) {
866 CHKERR moab.add_entities(meshset, *faces);
868 CHKERR moab.write_file(name.c_str(),
"VTK",
"", &meshset, 1, ths, 4);
869 CHKERR moab.delete_entities(&meshset, 1);
891 if (
B != PETSC_NULL) {
895 if (
F != PETSC_NULL) {
914 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
915 "Get edge sliding constrains element scaling",
918 CHKERR PetscOptionsScalar(
"-edge_sliding_alpha",
"scaling parameter",
"",
920 ierr = PetscOptionsEnd();
930 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
945 boost::shared_ptr<VectorDouble> &active_variables_ptr,
946 boost::shared_ptr<VectorDouble> &results_ptr,
947 boost::shared_ptr<MatrixDouble> &jacobian_ptr,
948 bool evaluate_jacobian,
const double &alpha)
958 if (
type != MBVERTEX)
966 Tag th0, th1, th2, th3;
980 ublas::vector<adouble> lambda_dofs(4);
981 for (
int dd = 0;
dd != 4; ++
dd) {
982 lambda_dofs[
dd] <<= (*activeVariablesPtr)[
dd];
984 ublas::vector<adouble> position_dofs(6);
985 for (
int dd = 0;
dd != 6; ++
dd) {
986 position_dofs[
dd] <<= (*activeVariablesPtr)[4 +
dd];
990 &position_dofs[0], &position_dofs[1], &position_dofs[2]);
992 &position_dofs[3], &position_dofs[4], &position_dofs[5]);
995 t_tangent(
i) = t_node1(
i) - t_node0(
i);
996 adouble l = sqrt(t_tangent(
i) * t_tangent(
i));
1000 t_dot0 = t_edge_base0(
i) * t_tangent(
i);
1001 t_dot1 = t_edge_base1(
i) * t_tangent(
i);
1004 t_base0(
i) = t_edge_base0(
i) - t_dot0 * t_tangent(
i);
1005 t_base1(
i) = t_edge_base1(
i) - t_dot1 * t_tangent(
i);
1006 t_base0(
i) /= sqrt(t_base0(
i) * t_base0(
i));
1007 t_base1(
i) /= sqrt(t_base1(
i) * t_base1(
i));
1009 auto t_base_fun1 = data.getFTensor0N();
1010 auto t_base_fun2 = data.getFTensor0N();
1016 ublas::vector<adouble> c_vec(4);
1017 ublas::vector<adouble> f_vec(6);
1021 int nb_gauss_pts = data.getN().size1();
1022 int nb_base_functions = data.getN().size2();
1023 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1026 &position_dofs[0], &position_dofs[1], &position_dofs[2], 3);
1028 &lambda_dofs[1], 2);
1032 for (
int bb = 0; bb != nb_base_functions; ++bb) {
1033 t_position(
i) += t_base_fun1 * t_position_dofs(
i);
1034 t_lambda(
j) += t_base_fun1 * t_lambda_dof(
j);
1040 t_delta(
i) = t_position(
i) - t_coord_ref(
i);
1048 for (
int bb = 0; bb != nb_base_functions; ++bb) {
1049 if (indices[2 * bb] != -1) {
1050 val =
w * t_base_fun2;
1051 t_c(N0) += val * dot0;
1052 t_c(N1) += val * dot1;
1053 val1 = val * t_lambda(N0);
1054 val2 = val * t_lambda(N1);
1055 t_f(
i) += val1 * t_base0(
i) + val2 * t_base1(
i);
1065 for (
int rr = 0; rr != 4; ++rr) {
1066 c_vec[rr] >>= (*resultsPtr)[rr];
1068 for (
int rr = 0; rr != 6; ++rr) {
1069 f_vec(rr) >>= (*resultsPtr)[4 + rr];
1075 double *jac_ptr[4 + 6];
1076 for (
int rr = 0; rr != 4 + 6; ++rr) {
1077 jac_ptr[rr] = &(*jacobianPtr)(rr, 0);
1084 "ADOL-C function evaluation with error");
1093 const std::string lagrange_multipliers_field_name,
1094 const std::string material_field_name) {
1099 material_field_name);
1104 const std::string lagrange_multipliers_field_name,
1105 const std::string material_field_name,
1106 const double *alpha =
nullptr) {
1113 boost::shared_ptr<VectorDouble> active_variables_ptr(
1115 boost::shared_ptr<VectorDouble> results_ptr(
new VectorDouble(4 + 6));
1116 boost::shared_ptr<MatrixDouble> jacobian_ptr(
1121 lagrange_multipliers_field_name, active_variables_ptr));
1123 material_field_name, active_variables_ptr));
1125 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1126 jacobian_ptr,
false,
aLpha));
1135 lagrange_multipliers_field_name, active_variables_ptr));
1137 material_field_name, active_variables_ptr));
1139 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1140 jacobian_ptr,
true,
aLpha));
1142 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1144 material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
1146 material_field_name, material_field_name, jacobian_ptr));
1153 const std::string lagrange_multipliers_field_name,
1154 const std::string material_field_name) {
1160 material_field_name);
1166 const std::string lagrange_multipliers_field_name,
1167 const std::string material_field_name) {
1170 boost::shared_ptr<VectorDouble> active_variables_ptr(
1172 boost::shared_ptr<VectorDouble> results_ptr(
new VectorDouble(4 + 6));
1173 boost::shared_ptr<MatrixDouble> jacobian_ptr(
1179 lagrange_multipliers_field_name, active_variables_ptr));
1181 material_field_name, active_variables_ptr));
1183 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1184 jacobian_ptr,
true,
aLpha));
1186 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1192 #endif // WITH_ADOL_C
1194 #endif // __SURFACE_SLIDING_CONSTRAINS_HPP__