29 boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET),
30 boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET),
31 boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET),
32 boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET),
33 boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET)
39 boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET),
40 boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET),
41 boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET),
42 boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET),
43 boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET)
49 boost::make_shared<DerivedDataForcesAndSourcesCore>(
51 boost::make_shared<DerivedDataForcesAndSourcesCore>(dataOnMaster[
H1]),
52 boost::make_shared<DerivedDataForcesAndSourcesCore>(
54 boost::make_shared<DerivedDataForcesAndSourcesCore>(
56 boost::make_shared<DerivedDataForcesAndSourcesCore>(dataOnMaster[
L2])
62 boost::make_shared<DerivedDataForcesAndSourcesCore>(
64 boost::make_shared<DerivedDataForcesAndSourcesCore>(dataOnSlave[
H1]),
65 boost::make_shared<DerivedDataForcesAndSourcesCore>(
67 boost::make_shared<DerivedDataForcesAndSourcesCore>(
69 boost::make_shared<DerivedDataForcesAndSourcesCore>(dataOnSlave[
L2])
72 dataH1Master(*dataOnMaster[
H1].get()),
73 dataH1Slave(*dataOnSlave[
H1].get()),
74 dataNoFieldSlave(*dataOnSlave[
NOFIELD].get()),
75 dataNoFieldMaster(*dataOnMaster[
NOFIELD].get()),
76 dataHcurlMaster(*dataOnMaster[
HCURL].get()),
77 dataHcurlSlave(*dataOnSlave[
HCURL].get()),
78 dataHdivMaster(*dataOnMaster[
HDIV].get()),
79 dataL2Master(*dataOnMaster[
L2].get()),
80 dataHdivSlave(*dataOnSlave[
HDIV].get()),
81 dataL2Slave(*dataOnSlave[
L2].get()), opContravariantTransform() {
83 getUserPolynomialBase() =
84 boost::shared_ptr<BaseFunction>(
new TriPolynomialBase());
86 dataH1Master.dataOnEntities[MBVERTEX].push_back(
88 dataH1Slave.dataOnEntities[MBVERTEX].push_back(
91 for (
int ee = 0; ee != 3; ++ee) {
92 dataH1Master.dataOnEntities[MBEDGE].push_back(
94 dataH1Slave.dataOnEntities[MBEDGE].push_back(
98 dataH1Master.dataOnEntities[MBTRI].push_back(
100 dataH1Slave.dataOnEntities[MBTRI].push_back(
103 dataHdivMaster.dataOnEntities[MBTRI].push_back(
105 dataHdivSlave.dataOnEntities[MBTRI].push_back(
109 dataOnMaster[
H1]->setElementType(MBTRI);
110 derivedDataOnMaster[
H1]->setElementType(MBTRI);
111 dataOnSlave[
H1]->setElementType(MBTRI);
112 derivedDataOnSlave[
H1]->setElementType(MBTRI);
114 dataOnMaster[
HDIV]->setElementType(MBTRI);
115 derivedDataOnMaster[
HDIV]->setElementType(MBTRI);
116 dataOnSlave[
HDIV]->setElementType(MBTRI);
117 derivedDataOnSlave[
HDIV]->setElementType(MBTRI);
152 double *shape_ptr_master =
155 shape_ptr_master, 1);
156 double *shape_ptr_slave =
181 auto get_coord_and_normal = [&]() {
186 coords.resize(num_nodes * 3,
false);
202 &vec_double(
r + 0), &vec_double(
r + 1), &vec_double(
r + 2));
205 auto t_coords_master = get_vec_ptr(
coords);
206 auto t_coords_slave = get_vec_ptr(
coords, 9);
207 auto t_normal_master = get_vec_ptr(
normal);
208 auto t_normal_slave = get_vec_ptr(
normal, 3);
217 &diff_ptr[0], &diff_ptr[1]);
225 for (
int nn = 0; nn != 3; ++nn) {
226 t_t1_master(
i) += t_coords_master(
i) * t_diff(
N0);
227 t_t1_slave(
i) += t_coords_slave(
i) * t_diff(
N0);
233 aRea[0] = sqrt(t_normal_master(
i) * t_normal_master(
i)) / 2.;
234 aRea[1] = sqrt(t_normal_slave(
i) * t_normal_slave(
i)) / 2.;
243 CHKERR get_coord_and_normal();
257 CHKERR getEntitySense<MBEDGE>(data_curl);
258 CHKERR getEntityDataOrder<MBEDGE>(data_curl,
HCURL);
259 CHKERR getEntitySense<MBTRI>(data_curl);
260 CHKERR getEntityDataOrder<MBTRI>(data_curl,
HCURL);
265 CHKERR getEntitySense<MBTRI>(data_div);
266 CHKERR getEntityDataOrder<MBTRI>(data_div,
HDIV);
274 CHKERR getEntitySense<MBTRI>(data_l2);
275 CHKERR getEntityDataOrder<MBTRI>(data_l2,
L2);
282 for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
283 data.spacesOnEntities[t].reset();
284 data.basesOnEntities[t].reset();
287 data.basesOnSpaces[s].reset();
296 if (shift != 0 && shift != 6) {
298 "Wrong shift for contact prism element");
301 data.
sPace = copy_data.sPace;
302 data.
bAse = copy_data.bAse;
303 for (
auto t : {MBVERTEX, MBEDGE, MBTRI}) {
309 for (
int ii = 0; ii != 3; ++ii) {
311 copy_data.dataOnEntities[MBEDGE][ii + shift].getSense();
313 copy_data.dataOnEntities[MBEDGE][ii + shift].getDataOrder();
318 copy_data.dataOnEntities[MBTRI][3].getSense();
320 copy_data.dataOnEntities[MBTRI][3].getDataOrder();
323 copy_data.dataOnEntities[MBTRI][4].getSense();
325 copy_data.dataOnEntities[MBTRI][4].getDataOrder();
336 if (shift != 3 && shift != 4) {
338 "Wrong shift for contact prism element");
341 data.
sPace = copy_data.sPace;
342 data.
bAse = copy_data.bAse;
344 for (
auto t : {MBVERTEX, MBTRI}) {
350 auto &cpy_ent_dat = copy_data.dataOnEntities[MBTRI][shift];
352 ent_dat.getBase() = cpy_ent_dat.getBase();
353 ent_dat.getSpace() = cpy_ent_dat.getSpace();
354 ent_dat.getSense() = ent_dat.getSense();
355 ent_dat.getDataOrder() = cpy_ent_dat.getDataOrder();
368 int rule =
getRule(order_row, order_col, order_data);
380 "Number of Gauss Points at Master triangle is different than slave");
397 CHKERR Tools::shapeFunMBTRI<1>(
411 for (
int dd = 0;
dd < 3;
dd++) {
473 &
m(0, 0), &
m(0, 1), &
m(0, 2));
502 "Not yet implemented");
506 "Not yet implemented");
511 "Not yet implemented");
526 constexpr std::array<UserDataOperator::OpType, 2> types{
528 std::array<std::string, 2> last_eval_field_name{std::string(), std::string()};
533 for (; oit != hi_oit; oit++) {
540 switch (oit->sPace) {
551 "Not implemented for this space", oit->sPace);
556 dataOnSlave[oit->sPace]->resetFieldDependentData();
566 boost::shared_ptr<DataForcesAndSourcesCore> op_master_data[2];
567 boost::shared_ptr<DataForcesAndSourcesCore> op_slave_data[2];
569 for (
int ss = 0; ss != 2; ss++) {
571 const std::string field_name =
572 !ss ? oit->rowFieldName : oit->colFieldName;
581 if ((oit->getNumeredEntFiniteElementPtr()->getBitFieldIdData() &
585 "no data field < %s > on finite element < %s >",
586 field_name.c_str(),
feName.c_str());
588 if (oit->getOpType() & types[ss] ||
603 "Not implemented for this space", space);
606 if (last_eval_field_name[ss] != field_name) {
612 boost::weak_ptr<EntityCacheNumeredDofs>
613 operator()(boost::shared_ptr<FieldEntity> &e) {
614 return e->entityCacheRowDofs;
620 MBPRISM, Extractor());
623 boost::weak_ptr<EntityCacheNumeredDofs>
624 operator()(boost::shared_ptr<FieldEntity> &e) {
625 return e->entityCacheColDofs;
630 MBPRISM, Extractor());
640 auto get_indices = [&](
auto &master,
auto &slave,
auto &ents,
644 master.dataOnEntities[MBVERTEX][0].getIndices(),
645 master.dataOnEntities[MBVERTEX][0].getLocalIndices(),
646 slave.dataOnEntities[MBVERTEX][0].getIndices(),
647 slave.dataOnEntities[MBVERTEX][0].getLocalIndices(), ex);
655 slave_data.dataOnEntities[MBVERTEX][0].getFieldData(),
657 slave_data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
659 slave_data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
661 slave_data.dataOnEntities[MBVERTEX][0].getSpace(),
663 slave_data.dataOnEntities[MBVERTEX][0].getBase());
669 boost::weak_ptr<EntityCacheNumeredDofs>
670 operator()(boost::shared_ptr<FieldEntity> &e) {
671 return e->entityCacheRowDofs;
675 CHKERR get_indices(*op_master_data[ss], *op_slave_data[ss],
679 boost::weak_ptr<EntityCacheNumeredDofs>
680 operator()(boost::shared_ptr<FieldEntity> &e) {
681 return e->entityCacheColDofs;
685 CHKERR get_indices(*op_master_data[ss], *op_slave_data[ss],
689 CHKERR get_data(*op_master_data[ss], *op_slave_data[ss]);
713 "Not implemented for this space", space);
715 last_eval_field_name[ss] = field_name;
724 type = cast_oit->getFaceType();
732 "Wrong combination of FaceType and OpType, OPROW or OPCOL "
733 "combined with face-face OpType");
740 "Wrong combination of FaceType and OpType, OPROWCOL "
741 "combined with face-face OpType");
746 "Face type is not set");
759 CHKERR oit->opRhs(*op_master_data[0],
false);
767 CHKERR oit->opRhs(*op_slave_data[0],
false);
775 CHKERR oit->opRhs(*op_master_data[1],
false);
783 CHKERR oit->opRhs(*op_slave_data[1],
false);
791 CHKERR oit->opLhs(*op_master_data[0], *op_master_data[1]);
799 CHKERR oit->opLhs(*op_master_data[0], *op_slave_data[1]);
807 CHKERR oit->opLhs(*op_slave_data[0], *op_master_data[1]);
815 CHKERR oit->opLhs(*op_slave_data[0], *op_slave_data[1]);
826 const std::string &field_name,
const EntityType type_lo,
827 const EntityType type_hi)
const {
830 auto reset_data = [type_lo, type_hi](
auto &data) {
831 for (EntityType t = type_lo; t != type_hi; ++t) {
832 for (
auto &dat : data.dataOnEntities[t]) {
833 dat.getDataOrder() = 0;
836 dat.getFieldData().resize(0,
false);
837 dat.getFieldDofs().resize(0,
false);
838 dat.getFieldEntities().resize(0,
false);
842 reset_data(master_data);
843 reset_data(slave_data);
849 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
851 if (lo != field_ents.end()) {
854 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid,
cmp_uid_hi);
856 for (
auto it = lo; it != hi; ++it)
857 if (
auto e = it->lock()) {
859 auto get_data = [&](
auto &data,
auto type,
auto side) {
860 auto &dat = data.dataOnEntities[
type][side];
861 auto &ent_field_dofs = dat.getFieldDofs();
862 auto &ent_field_data = dat.getFieldData();
863 dat.getFieldEntities().resize(1,
false);
864 dat.getFieldEntities()[0] = e.get();
865 dat.getBase() = e->getApproxBase();
866 dat.getSpace() = e->getSpace();
867 const int ent_order = e->getMaxOrder();
869 dat.getDataOrder() > ent_order ? dat.getDataOrder() : ent_order;
870 const auto dof_ent_field_data = e->getEntFieldData();
871 const int nb_dofs_on_ent = e->getNbDofsOnEnt();
872 ent_field_data.resize(nb_dofs_on_ent,
false);
873 noalias(ent_field_data) = e->getEntFieldData();
874 ent_field_dofs.resize(nb_dofs_on_ent,
false);
875 std::fill(ent_field_dofs.begin(), ent_field_dofs.end(),
nullptr);
876 if (
auto cache = e->entityCacheDataDofs.lock()) {
877 for (
auto dit =
cache->loHi[0]; dit !=
cache->loHi[1]; ++dit) {
878 ent_field_dofs[(*dit)->getEntDofIdx()] =
884 const EntityType
type = e->getEntType();
885 const int side = e->getSideNumberPtr()->side_number;
891 get_data(master_data,
type, side);
893 get_data(slave_data,
type, side - 6);
899 get_data(master_data,
type, 0);
901 get_data(slave_data,
type, 0);
906 "Entity type not implemented (FIXME)");
909 const int brother_side = e->getSideNumberPtr()->brother_side_number;
910 if (brother_side != -1)
912 "Case with brother side not implemented (FIXME)");
921 const std::string field_name,
VectorDouble &master_nodes_data,
932 auto set_zero = [](
auto &nodes_data,
auto &nodes_dofs,
auto &field_entities) {
933 nodes_data.resize(0,
false);
934 nodes_dofs.resize(0,
false);
935 field_entities.resize(0,
false);
937 set_zero(master_nodes_data, master_nodes_dofs, master_field_entities);
938 set_zero(slave_nodes_data, slave_nodes_dofs, slave_field_entities);
943 auto bit_number = (*field_it)->getBitNumber();
944 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
945 master_space = slave_space = (*field_it)->getSpace();
946 master_base = slave_base = (*field_it)->getApproxBase();
950 bit_number, get_id_for_min_type<MBVERTEX>());
951 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
953 if (lo != field_ents.end()) {
955 bit_number, get_id_for_max_type<MBVERTEX>());
956 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid,
cmp_uid_hi);
960 for (
auto it = lo; it != hi; ++it) {
961 if (
auto e = it->lock()) {
962 nb_dofs += e->getEntFieldData().size();
968 auto init_set = [&](
auto &nodes_data,
auto &nodes_dofs,
969 auto &field_entities) {
970 constexpr
int num_nodes = 3;
971 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
972 nodes_data.resize(max_nb_dofs,
false);
973 nodes_dofs.resize(max_nb_dofs,
false);
974 field_entities.resize(num_nodes,
false);
976 fill(nodes_dofs.begin(), nodes_dofs.end(),
nullptr);
977 fill(field_entities.begin(), field_entities.end(),
nullptr);
980 init_set(master_nodes_data, master_nodes_dofs, master_field_entities);
981 init_set(slave_nodes_data, slave_nodes_dofs, slave_field_entities);
983 for (
auto it = lo; it != hi; ++it) {
984 if (
auto e = it->lock()) {
986 const auto &sn = e->getSideNumberPtr();
987 int side = sn->side_number;
989 auto set_data = [&](
auto &nodes_data,
auto &nodes_dofs,
990 auto &field_entities,
int side,
int pos) {
991 field_entities[side] = e.get();
992 if (
auto cache = e->entityCacheDataDofs.lock()) {
993 for (
auto dit =
cache->loHi[0]; dit !=
cache->loHi[1];
995 const auto dof_idx = (*dit)->getEntDofIdx();
996 nodes_data[pos + dof_idx] = (*dit)->getFieldData();
997 nodes_dofs[pos + dof_idx] =
1004 set_data(master_nodes_data, master_nodes_dofs,
1005 master_field_entities, side, side * nb_dofs_on_vert);
1007 set_data(slave_nodes_data, slave_nodes_dofs,
1008 slave_field_entities, (side - 3),
1009 (side - 3) * nb_dofs_on_vert);
1011 const int brother_side = sn->brother_side_number;
1012 if (brother_side != -1)
1014 "Not implemented (FIXME please)");
1025 template <
typename EXTRACTOR>
1029 const EntityType type_lo,
const EntityType type_hi,
1030 EXTRACTOR &&extractor)
const {
1033 auto clear_data = [type_lo, type_hi](
auto &data) {
1034 for (EntityType t = type_lo; t != type_hi; ++t) {
1035 for (
auto &
d : data.dataOnEntities[t]) {
1036 d.getIndices().resize(0,
false);
1037 d.getLocalIndices().resize(0,
false);
1042 clear_data(master_data);
1043 clear_data(slave_data);
1048 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
1050 if (lo != ents_field.end()) {
1053 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid,
cmp_uid_hi);
1056 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1058 for (
auto it = lo; it != hi; ++it)
1059 if (
auto e = it->lock()) {
1061 const EntityType
type = e->getEntType();
1062 const int side = e->getSideNumberPtr()->side_number;
1063 const int brother_side = e->getSideNumberPtr()->brother_side_number;
1064 if (brother_side != -1)
1066 "Not implemented case");
1068 auto get_indices = [&](
auto &data,
const auto type,
const auto side) {
1069 if (
auto cache = extractor(e).lock()) {
1070 for (
auto dit =
cache->loHi[0]; dit !=
cache->loHi[1]; ++dit) {
1071 auto &dof = (**dit);
1072 auto &dat = data.dataOnEntities[
type][side];
1073 auto &ent_field_indices = dat.getIndices();
1074 auto &ent_field_local_indices = dat.getLocalIndices();
1075 if (ent_field_indices.empty()) {
1076 const int nb_dofs_on_ent = dof.getNbDofsOnEnt();
1077 ent_field_indices.resize(nb_dofs_on_ent,
false);
1078 ent_field_local_indices.resize(nb_dofs_on_ent,
false);
1079 std::fill(ent_field_indices.data().begin(),
1080 ent_field_indices.data().end(), -1);
1081 std::fill(ent_field_local_indices.data().begin(),
1082 ent_field_local_indices.data().end(), -1);
1084 const int idx = dof.getEntDofIdx();
1085 ent_field_indices[idx] = dof.getPetscGlobalDofIdx();
1086 ent_field_local_indices[idx] = dof.getPetscLocalDofIdx();
1095 get_indices(master_data,
type, side);
1097 get_indices(slave_data,
type, side - 6);
1103 get_indices(master_data,
type, 0);
1105 get_indices(slave_data,
type, 0);
1110 "Entity type not implemented");
1119 template <
typename EXTRACTOR>
1124 EXTRACTOR &&extractor)
const {
1127 master_nodes_indices.resize(0,
false);
1128 master_local_nodes_indices.resize(0,
false);
1129 slave_nodes_indices.resize(0,
false);
1130 slave_local_nodes_indices.resize(0,
false);
1135 auto bit_number = (*field_it)->getBitNumber();
1136 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
1139 bit_number, get_id_for_min_type<MBVERTEX>());
1140 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
1142 if (lo != ents_field.end()) {
1144 bit_number, get_id_for_max_type<MBVERTEX>());
1145 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid,
cmp_uid_hi);
1149 for (
auto it = lo; it != hi; ++it) {
1150 if (
auto e = it->lock()) {
1151 if (
auto cache = extractor(e).lock()) {
1152 nb_dofs += std::distance(
cache->loHi[0],
cache->loHi[1]);
1159 constexpr
int num_nodes = 3;
1160 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
1162 auto set_vec_size = [&](
auto &nodes_indices,
1163 auto &local_nodes_indices) {
1164 nodes_indices.resize(max_nb_dofs,
false);
1165 local_nodes_indices.resize(max_nb_dofs,
false);
1166 std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
1167 std::fill(local_nodes_indices.begin(), local_nodes_indices.end(),
1171 set_vec_size(master_nodes_indices, master_local_nodes_indices);
1172 set_vec_size(slave_nodes_indices, slave_local_nodes_indices);
1174 for (
auto it = lo; it != hi; ++it) {
1175 if (
auto e = it->lock()) {
1177 const int side = e->getSideNumberPtr()->side_number;
1178 const int brother_side =
1179 e->getSideNumberPtr()->brother_side_number;
1180 if (brother_side != -1)
1182 "Not implemented case");
1184 auto get_indices = [&](
auto &nodes_indices,
1185 auto &local_nodes_indices,
1187 if (
auto cache = extractor(e).lock()) {
1188 for (
auto dit =
cache->loHi[0]; dit !=
cache->loHi[1];
1190 const int idx = (*dit)->getPetscGlobalDofIdx();
1191 const int local_idx = (*dit)->getPetscLocalDofIdx();
1193 side * nb_dofs_on_vert + (*dit)->getDofCoeffIdx();
1194 nodes_indices[pos] = idx;
1195 local_nodes_indices[pos] = local_idx;
1201 get_indices(master_nodes_indices, master_local_nodes_indices,
1204 get_indices(slave_nodes_indices, slave_local_nodes_indices,