17 boost::make_shared<EntitiesFieldData>(MBENTITYSET),
18 boost::make_shared<EntitiesFieldData>(MBENTITYSET),
19 boost::make_shared<EntitiesFieldData>(MBENTITYSET),
20 boost::make_shared<EntitiesFieldData>(MBENTITYSET),
21 boost::make_shared<EntitiesFieldData>(MBENTITYSET)
27 boost::make_shared<EntitiesFieldData>(MBENTITYSET),
28 boost::make_shared<EntitiesFieldData>(MBENTITYSET),
29 boost::make_shared<EntitiesFieldData>(MBENTITYSET),
30 boost::make_shared<EntitiesFieldData>(MBENTITYSET),
31 boost::make_shared<EntitiesFieldData>(MBENTITYSET)
37 boost::make_shared<DerivedEntitiesFieldData>(
39 boost::make_shared<DerivedEntitiesFieldData>(dataOnMaster[
H1]),
40 boost::make_shared<DerivedEntitiesFieldData>(
42 boost::make_shared<DerivedEntitiesFieldData>(
44 boost::make_shared<DerivedEntitiesFieldData>(dataOnMaster[
L2])
50 boost::make_shared<DerivedEntitiesFieldData>(
52 boost::make_shared<DerivedEntitiesFieldData>(dataOnSlave[
H1]),
53 boost::make_shared<DerivedEntitiesFieldData>(
55 boost::make_shared<DerivedEntitiesFieldData>(
57 boost::make_shared<DerivedEntitiesFieldData>(dataOnSlave[
L2])
60 dataH1Master(*dataOnMaster[
H1].get()),
61 dataH1Slave(*dataOnSlave[
H1].get()),
62 dataNoFieldSlave(*dataOnSlave[
NOFIELD].get()),
63 dataNoFieldMaster(*dataOnMaster[
NOFIELD].get()),
64 dataHcurlMaster(*dataOnMaster[
HCURL].get()),
65 dataHcurlSlave(*dataOnSlave[
HCURL].get()),
66 dataHdivMaster(*dataOnMaster[
HDIV].get()),
67 dataL2Master(*dataOnMaster[
L2].get()),
68 dataHdivSlave(*dataOnSlave[
HDIV].get()),
69 dataL2Slave(*dataOnSlave[
L2].get()), opContravariantTransform() {
71 getUserPolynomialBase() =
72 boost::shared_ptr<BaseFunction>(
new TriPolynomialBase());
75 dataOnMaster[
H1]->setElementType(MBTRI);
76 derivedDataOnMaster[
H1]->setElementType(MBTRI);
77 dataOnSlave[
H1]->setElementType(MBTRI);
78 derivedDataOnSlave[
H1]->setElementType(MBTRI);
80 dataOnMaster[
HDIV]->setElementType(MBTRI);
81 derivedDataOnMaster[
HDIV]->setElementType(MBTRI);
82 dataOnSlave[
HDIV]->setElementType(MBTRI);
83 derivedDataOnSlave[
HDIV]->setElementType(MBTRI);
85 dataOnMaster[
NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
87 dataOnSlave[
NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
90 derivedDataOnMaster[
NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
92 derivedDataOnSlave[
NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
96 "Problem with creation data on element");
132 double *shape_ptr_master =
135 shape_ptr_master, 1);
136 double *shape_ptr_slave =
160 auto get_coord_and_normal = [&]() {
165 coords.resize(num_nodes * 3,
false);
181 &vec_double(
r + 0), &vec_double(
r + 1), &vec_double(
r + 2));
184 auto t_coords_master = get_vec_ptr(
coords);
185 auto t_coords_slave = get_vec_ptr(
coords, 9);
186 auto t_normal_master = get_vec_ptr(
normal);
187 auto t_normal_slave = get_vec_ptr(
normal, 3);
196 &diff_ptr[0], &diff_ptr[1]);
204 for (
int nn = 0; nn != 3; ++nn) {
205 t_t1_master(
i) += t_coords_master(
i) * t_diff(N0);
206 t_t1_slave(
i) += t_coords_slave(
i) * t_diff(N0);
212 aRea[0] = sqrt(t_normal_master(
i) * t_normal_master(
i)) / 2.;
213 aRea[1] = sqrt(t_normal_slave(
i) * t_normal_slave(
i)) / 2.;
222 CHKERR get_coord_and_normal();
236 CHKERR getEntitySense<MBEDGE>(data_curl);
237 CHKERR getEntityDataOrder<MBEDGE>(data_curl,
HCURL);
238 CHKERR getEntitySense<MBTRI>(data_curl);
239 CHKERR getEntityDataOrder<MBTRI>(data_curl,
HCURL);
244 CHKERR getEntitySense<MBTRI>(data_div);
245 CHKERR getEntityDataOrder<MBTRI>(data_div,
HDIV);
253 CHKERR getEntitySense<MBTRI>(data_l2);
254 CHKERR getEntityDataOrder<MBTRI>(data_l2,
L2);
261 for (EntityType
t = MBVERTEX;
t != MBMAXTYPE; ++
t) {
262 data.spacesOnEntities[
t].reset();
263 data.basesOnEntities[
t].reset();
266 data.basesOnSpaces[s].reset();
275 if (shift != 0 && shift != 6) {
277 "Wrong shift for contact prism element");
280 data.
sPace = copy_data.sPace;
281 data.
bAse = copy_data.bAse;
282 for (
auto t : {MBVERTEX, MBEDGE, MBTRI}) {
288 for (
int ii = 0; ii != 3; ++ii) {
290 copy_data.dataOnEntities[MBEDGE][ii + shift].getSense();
292 copy_data.dataOnEntities[MBEDGE][ii + shift].getOrder();
297 copy_data.dataOnEntities[MBTRI][3].getSense();
299 copy_data.dataOnEntities[MBTRI][3].getOrder();
302 copy_data.dataOnEntities[MBTRI][4].getSense();
304 copy_data.dataOnEntities[MBTRI][4].getOrder();
315 if (shift != 3 && shift != 4) {
317 "Wrong shift for contact prism element");
320 data.
sPace = copy_data.sPace;
321 data.
bAse = copy_data.bAse;
323 for (
auto t : {MBVERTEX, MBTRI}) {
329 auto &cpy_ent_dat = copy_data.dataOnEntities[MBTRI][shift];
331 ent_dat.getBase() = cpy_ent_dat.getBase();
332 ent_dat.getSpace() = cpy_ent_dat.getSpace();
333 ent_dat.getSense() = ent_dat.getSense();
334 ent_dat.getOrder() = cpy_ent_dat.getOrder();
347 int rule =
getRule(order_row, order_col, order_data);
359 "Number of Gauss Points at Master triangle is different than slave");
376 CHKERR Tools::shapeFunMBTRI<1>(
390 for (
int dd = 0;
dd < 3;
dd++) {
471 "Not yet implemented");
475 "Not yet implemented");
480 "Not yet implemented");
494 constexpr std::array<UserDataOperator::OpType, 2> types{
496 std::array<std::string, 2> last_eval_field_name{std::string(), std::string()};
501 for (; oit != hi_oit; oit++) {
508 switch (oit->sPace) {
519 "Not implemented for this space", oit->sPace);
524 dataOnSlave[oit->sPace]->resetFieldDependentData();
534 boost::shared_ptr<EntitiesFieldData> op_master_data[2];
535 boost::shared_ptr<EntitiesFieldData> op_slave_data[2];
537 for (
int ss = 0; ss != 2; ss++) {
540 !ss ? oit->rowFieldName : oit->colFieldName;
549 if ((oit->getNumeredEntFiniteElementPtr()->getBitFieldIdData() &
553 "no data field < %s > on finite element < %s >",
556 if (oit->getOpType() & types[ss] ||
571 "Not implemented for this space", space);
580 boost::weak_ptr<EntityCacheNumeredDofs>
581 operator()(boost::shared_ptr<FieldEntity> &e) {
582 return e->entityCacheRowDofs;
588 MBPRISM, Extractor());
591 boost::weak_ptr<EntityCacheNumeredDofs>
592 operator()(boost::shared_ptr<FieldEntity> &e) {
593 return e->entityCacheColDofs;
598 MBPRISM, Extractor());
608 auto get_indices = [&](
auto &master,
auto &slave,
auto &ents,
612 master.dataOnEntities[MBVERTEX][0].getIndices(),
613 master.dataOnEntities[MBVERTEX][0].getLocalIndices(),
614 slave.dataOnEntities[MBVERTEX][0].getIndices(),
615 slave.dataOnEntities[MBVERTEX][0].getLocalIndices(), ex);
623 slave_data.dataOnEntities[MBVERTEX][0].getFieldData(),
625 slave_data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
627 slave_data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
629 slave_data.dataOnEntities[MBVERTEX][0].getSpace(),
631 slave_data.dataOnEntities[MBVERTEX][0].getBase());
637 boost::weak_ptr<EntityCacheNumeredDofs>
638 operator()(boost::shared_ptr<FieldEntity> &e) {
639 return e->entityCacheRowDofs;
643 CHKERR get_indices(*op_master_data[ss], *op_slave_data[ss],
647 boost::weak_ptr<EntityCacheNumeredDofs>
648 operator()(boost::shared_ptr<FieldEntity> &e) {
649 return e->entityCacheColDofs;
653 CHKERR get_indices(*op_master_data[ss], *op_slave_data[ss],
657 CHKERR get_data(*op_master_data[ss], *op_slave_data[ss]);
682 "Not implemented for this space", space);
693 type = cast_oit->getFaceType();
701 "Wrong combination of FaceType and OpType, OPROW or OPCOL "
702 "combined with face-face OpType");
709 "Wrong combination of FaceType and OpType, OPROWCOL "
710 "combined with face-face OpType");
715 "Face type is not set");
728 CHKERR oit->opRhs(*op_master_data[0],
false);
736 CHKERR oit->opRhs(*op_slave_data[0],
false);
744 CHKERR oit->opRhs(*op_master_data[1],
false);
752 CHKERR oit->opRhs(*op_slave_data[1],
false);
760 CHKERR oit->opLhs(*op_master_data[0], *op_master_data[1]);
768 CHKERR oit->opLhs(*op_master_data[0], *op_slave_data[1]);
776 CHKERR oit->opLhs(*op_slave_data[0], *op_master_data[1]);
784 CHKERR oit->opLhs(*op_slave_data[0], *op_slave_data[1]);
795 const std::string &
field_name,
const EntityType type_lo,
796 const EntityType type_hi)
const {
799 auto reset_data = [type_lo, type_hi](
auto &data) {
800 for (EntityType
t = type_lo;
t != type_hi; ++
t) {
801 for (
auto &dat : data.dataOnEntities[
t]) {
805 dat.getFieldData().resize(0,
false);
806 dat.getFieldDofs().resize(0,
false);
807 dat.getFieldEntities().resize(0,
false);
811 reset_data(master_data);
812 reset_data(slave_data);
818 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
820 if (lo != field_ents.end()) {
823 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid,
cmp_uid_hi);
825 for (
auto it = lo; it != hi; ++it)
826 if (
auto e = it->lock()) {
828 auto get_data = [&](
auto &data,
auto type,
auto side) {
829 auto &dat = data.dataOnEntities[
type][side];
830 auto &ent_field_dofs = dat.getFieldDofs();
831 auto &ent_field_data = dat.getFieldData();
832 dat.getFieldEntities().resize(1,
false);
833 dat.getFieldEntities()[0] = e.get();
834 dat.getBase() = e->getApproxBase();
835 dat.getSpace() = e->getSpace();
836 const int ent_order = e->getMaxOrder();
838 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
839 const auto dof_ent_field_data = e->getEntFieldData();
840 const int nb_dofs_on_ent = e->getNbDofsOnEnt();
841 ent_field_data.resize(nb_dofs_on_ent,
false);
842 noalias(ent_field_data) = e->getEntFieldData();
843 ent_field_dofs.resize(nb_dofs_on_ent,
false);
844 std::fill(ent_field_dofs.begin(), ent_field_dofs.end(),
nullptr);
845 if (
auto cache = e->entityCacheDataDofs.lock()) {
846 for (
auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
847 ent_field_dofs[(*dit)->getEntDofIdx()] =
853 const EntityType
type = e->getEntType();
854 const int side = e->getSideNumberPtr()->side_number;
860 get_data(master_data,
type, side);
862 get_data(slave_data,
type, side - 6);
868 get_data(master_data,
type, 0);
870 get_data(slave_data,
type, 0);
875 "Entity type not implemented (FIXME)");
878 const int brother_side = e->getSideNumberPtr()->brother_side_number;
879 if (brother_side != -1)
881 "Case with brother side not implemented (FIXME)");
901 auto set_zero = [](
auto &nodes_data,
auto &nodes_dofs,
auto &field_entities) {
902 nodes_data.resize(0,
false);
903 nodes_dofs.resize(0,
false);
904 field_entities.resize(0,
false);
906 set_zero(master_nodes_data, master_nodes_dofs, master_field_entities);
907 set_zero(slave_nodes_data, slave_nodes_dofs, slave_field_entities);
912 auto bit_number = (*field_it)->getBitNumber();
913 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
914 master_space = slave_space = (*field_it)->getSpace();
915 master_base = slave_base = (*field_it)->getApproxBase();
919 bit_number, get_id_for_min_type<MBVERTEX>());
920 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
922 if (lo != field_ents.end()) {
924 bit_number, get_id_for_max_type<MBVERTEX>());
925 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid,
cmp_uid_hi);
929 for (
auto it = lo; it != hi; ++it) {
930 if (
auto e = it->lock()) {
931 nb_dofs += e->getEntFieldData().size();
937 auto init_set = [&](
auto &nodes_data,
auto &nodes_dofs,
938 auto &field_entities) {
939 constexpr
int num_nodes = 3;
940 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
941 nodes_data.resize(max_nb_dofs,
false);
942 nodes_dofs.resize(max_nb_dofs,
false);
943 field_entities.resize(num_nodes,
false);
945 fill(nodes_dofs.begin(), nodes_dofs.end(),
nullptr);
946 fill(field_entities.begin(), field_entities.end(),
nullptr);
949 init_set(master_nodes_data, master_nodes_dofs, master_field_entities);
950 init_set(slave_nodes_data, slave_nodes_dofs, slave_field_entities);
952 for (
auto it = lo; it != hi; ++it) {
953 if (
auto e = it->lock()) {
955 const auto &sn = e->getSideNumberPtr();
956 int side = sn->side_number;
958 auto set_data = [&](
auto &nodes_data,
auto &nodes_dofs,
959 auto &field_entities,
int side,
int pos) {
960 field_entities[side] = e.get();
961 if (
auto cache = e->entityCacheDataDofs.lock()) {
962 for (
auto dit = cache->loHi[0]; dit != cache->loHi[1];
964 const auto dof_idx = (*dit)->getEntDofIdx();
965 nodes_data[pos + dof_idx] = (*dit)->getFieldData();
966 nodes_dofs[pos + dof_idx] =
973 set_data(master_nodes_data, master_nodes_dofs,
974 master_field_entities, side, side * nb_dofs_on_vert);
976 set_data(slave_nodes_data, slave_nodes_dofs,
977 slave_field_entities, (side - 3),
978 (side - 3) * nb_dofs_on_vert);
980 const int brother_side = sn->brother_side_number;
981 if (brother_side != -1)
983 "Not implemented (FIXME please)");
994 template <
typename EXTRACTOR>
998 const EntityType type_lo,
const EntityType type_hi,
999 EXTRACTOR &&extractor)
const {
1002 auto clear_data = [type_lo, type_hi](
auto &data) {
1003 for (EntityType
t = type_lo;
t != type_hi; ++
t) {
1004 for (
auto &
d : data.dataOnEntities[
t]) {
1005 d.getIndices().resize(0,
false);
1006 d.getLocalIndices().resize(0,
false);
1011 clear_data(master_data);
1012 clear_data(slave_data);
1017 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
1019 if (lo != ents_field.end()) {
1022 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid,
cmp_uid_hi);
1025 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1027 for (
auto it = lo; it != hi; ++it)
1028 if (
auto e = it->lock()) {
1030 const EntityType
type = e->getEntType();
1031 const int side = e->getSideNumberPtr()->side_number;
1032 const int brother_side = e->getSideNumberPtr()->brother_side_number;
1033 if (brother_side != -1)
1035 "Not implemented case");
1037 auto get_indices = [&](
auto &data,
const auto type,
const auto side) {
1038 if (
auto cache = extractor(e).lock()) {
1039 for (
auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
1040 auto &dof = (**dit);
1041 auto &dat = data.dataOnEntities[
type][side];
1042 auto &ent_field_indices = dat.getIndices();
1043 auto &ent_field_local_indices = dat.getLocalIndices();
1044 if (ent_field_indices.empty()) {
1045 const int nb_dofs_on_ent = dof.getNbDofsOnEnt();
1046 ent_field_indices.resize(nb_dofs_on_ent,
false);
1047 ent_field_local_indices.resize(nb_dofs_on_ent,
false);
1048 std::fill(ent_field_indices.data().begin(),
1049 ent_field_indices.data().end(), -1);
1050 std::fill(ent_field_local_indices.data().begin(),
1051 ent_field_local_indices.data().end(), -1);
1053 const int idx = dof.getEntDofIdx();
1054 ent_field_indices[idx] = dof.getPetscGlobalDofIdx();
1055 ent_field_local_indices[idx] = dof.getPetscLocalDofIdx();
1064 get_indices(master_data,
type, side);
1066 get_indices(slave_data,
type, side - 6);
1072 get_indices(master_data,
type, 0);
1074 get_indices(slave_data,
type, 0);
1079 "Entity type not implemented");
1088 template <
typename EXTRACTOR>
1093 EXTRACTOR &&extractor)
const {
1096 master_nodes_indices.resize(0,
false);
1097 master_local_nodes_indices.resize(0,
false);
1098 slave_nodes_indices.resize(0,
false);
1099 slave_local_nodes_indices.resize(0,
false);
1104 auto bit_number = (*field_it)->getBitNumber();
1105 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
1108 bit_number, get_id_for_min_type<MBVERTEX>());
1109 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
1111 if (lo != ents_field.end()) {
1113 bit_number, get_id_for_max_type<MBVERTEX>());
1114 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid,
cmp_uid_hi);
1118 for (
auto it = lo; it != hi; ++it) {
1119 if (
auto e = it->lock()) {
1120 if (
auto cache = extractor(e).lock()) {
1121 nb_dofs += std::distance(cache->loHi[0], cache->loHi[1]);
1128 constexpr
int num_nodes = 3;
1129 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
1131 auto set_vec_size = [&](
auto &nodes_indices,
1132 auto &local_nodes_indices) {
1133 nodes_indices.resize(max_nb_dofs,
false);
1134 local_nodes_indices.resize(max_nb_dofs,
false);
1135 std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
1136 std::fill(local_nodes_indices.begin(), local_nodes_indices.end(),
1140 set_vec_size(master_nodes_indices, master_local_nodes_indices);
1141 set_vec_size(slave_nodes_indices, slave_local_nodes_indices);
1143 for (
auto it = lo; it != hi; ++it) {
1144 if (
auto e = it->lock()) {
1146 const int side = e->getSideNumberPtr()->side_number;
1147 const int brother_side =
1148 e->getSideNumberPtr()->brother_side_number;
1149 if (brother_side != -1)
1151 "Not implemented case");
1153 auto get_indices = [&](
auto &nodes_indices,
1154 auto &local_nodes_indices,
1156 if (
auto cache = extractor(e).lock()) {
1157 for (
auto dit = cache->loHi[0]; dit != cache->loHi[1];
1159 const int idx = (*dit)->getPetscGlobalDofIdx();
1160 const int local_idx = (*dit)->getPetscLocalDofIdx();
1162 side * nb_dofs_on_vert + (*dit)->getDofCoeffIdx();
1163 nodes_indices[pos] = idx;
1164 local_nodes_indices[pos] = local_idx;
1170 get_indices(master_nodes_indices, master_local_nodes_indices,
1173 get_indices(slave_nodes_indices, slave_local_nodes_indices,
1190 const string fe_name,
1192 const int side_type,
const EntityHandle ent_for_side) {
1193 return loopSide(fe_name, &fe_method, side_type, ent_for_side);