33 {
37 int nb_gauss_pts = pts.size2();
38
40 auto &ptr = data.getBBAlphaIndicesSharedPtr(
field_name);
41 if (!ptr)
43 return *ptr;
44 };
45
48 if (!ptr)
50 return *ptr;
51 };
52
54 auto &ptr = data.getBBDiffNSharedPtr(
field_name);
55 if (!ptr)
57 return *ptr;
58 };
59
60 auto get_alpha_by_name_ptr =
61 [](auto &data,
62 const std::string &
field_name) -> boost::shared_ptr<MatrixInt> & {
63 return data.getBBAlphaIndicesSharedPtr(
field_name);
64 };
65
66 auto get_base_by_name_ptr =
67 [](auto &data,
68 const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
70 };
71
72 auto get_diff_base_by_name_ptr =
73 [](auto &data,
74 const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
76 };
77
78 auto get_alpha_by_order_ptr =
79 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
80 return data.getBBAlphaIndicesByOrderSharedPtr(o);
81 };
82
83 auto get_base_by_order_ptr =
84 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
85 return data.getBBNByOrderSharedPtr(o);
86 };
87
88 auto get_diff_base_by_order_ptr =
89 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
90 return data.getBBDiffNByOrderSharedPtr(o);
91 };
92
94 auto &vert_get_diff_n = get_diff_base(data.
dataOnEntities[MBVERTEX][0]);
95 vert_get_n.resize(nb_gauss_pts, 3, false);
96 vert_get_diff_n.resize(nb_gauss_pts, 6, false);
97
99 (unsigned int)nb_gauss_pts)
101 "Base functions or nodes has wrong number of integration points "
102 "for base %s",
105
107 vertex_alpha.resize(3, 3, false);
108 vertex_alpha.clear();
109 for (
int n = 0;
n != 3; ++
n)
111
113 1,
lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
115 &vert_get_diff_n(0, 0));
116 for (
int n = 0;
n != 3; ++
n) {
117 const int f = boost::math::factorial<double>(
119 for (
int g = 0;
g != nb_gauss_pts; ++
g) {
120 vert_get_n(
g,
n) *=
f;
121 for (
int d = 0;
d != 2; ++
d)
122 vert_get_diff_n(
g, 2 *
n + d) *=
f;
123 }
124 }
125
126
130 "Wrong size of ent data");
131
132 constexpr int edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
133 for (int ee = 0; ee != 3; ++ee) {
135
136 if (ent_data.getSense() == 0)
138 "Sense of the edge unknown");
139
140 const int sense = ent_data.getSense();
141 const int order = ent_data.getOrder();
142 const int nb_dofs =
NBEDGE_H1(ent_data.getOrder());
143
144 if (nb_dofs) {
145 if (get_alpha_by_order_ptr(ent_data,
order)) {
147 get_alpha_by_order_ptr(ent_data,
order);
149 get_base_by_order_ptr(ent_data,
order);
150 get_diff_base_by_name_ptr(ent_data,
field_name) =
151 get_diff_base_by_order_ptr(ent_data,
order);
152 } else {
153 auto &get_n = get_base(ent_data);
154 auto &get_diff_n = get_diff_base(ent_data);
155 get_n.resize(nb_gauss_pts, nb_dofs, false);
156 get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
157
159 edge_alpha.resize(nb_dofs, 3, false);
161 &edge_alpha(0, 0));
162 if (sense == -1) {
163 for (
int i = 0;
i != edge_alpha.size1(); ++
i) {
164 int a = edge_alpha(
i, edges_nodes[ee][0]);
165 edge_alpha(
i, edges_nodes[ee][0]) =
166 edge_alpha(
i, edges_nodes[ee][1]);
167 edge_alpha(
i, edges_nodes[ee][1]) =
a;
168 }
169 }
171 order,
lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
173 &get_diff_n(0, 0));
174
175 get_alpha_by_order_ptr(ent_data,
order) =
177 get_base_by_order_ptr(ent_data,
order) =
179 get_diff_base_by_order_ptr(ent_data,
order) =
180 get_diff_base_by_name_ptr(ent_data,
field_name);
181 }
182 }
183 }
184 } else {
185 for (int ee = 0; ee != 3; ++ee) {
187 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
188 auto &get_n = get_base(ent_data);
189 auto &get_diff_n = get_diff_base(ent_data);
190 get_n.resize(nb_gauss_pts, 0, false);
191 get_diff_n.resize(nb_gauss_pts, 0, false);
192 }
193 }
194
198 "Wrong size ent of ent data");
199
201 const int order = ent_data.getOrder();
203 if (get_alpha_by_order_ptr(ent_data,
order)) {
205 get_alpha_by_order_ptr(ent_data,
order);
207 get_base_by_order_ptr(ent_data,
order);
208 get_diff_base_by_name_ptr(ent_data,
field_name) =
209 get_diff_base_by_order_ptr(ent_data,
order);
210 } else {
211
212 auto &get_n = get_base(ent_data);
213 auto &get_diff_n = get_diff_base(ent_data);
214 get_n.resize(nb_gauss_pts, nb_dofs, false);
215 get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
216 if (nb_dofs) {
217 auto &tri_alpha = get_alpha(ent_data);
218 tri_alpha.resize(nb_dofs, 3, false);
219
222 order,
lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
224 &get_diff_n(0, 0));
225
226 get_alpha_by_order_ptr(ent_data,
order) =
228 get_base_by_order_ptr(ent_data,
order) =
230 get_diff_base_by_order_ptr(ent_data,
order) =
231 get_diff_base_by_name_ptr(ent_data,
field_name);
232 }
233 }
234 } else {
236 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
237 auto &get_n = get_base(ent_data);
238 auto &get_diff_n = get_diff_base(ent_data);
239 get_n.resize(nb_gauss_pts, 0, false);
240 get_diff_n.resize(nb_gauss_pts, 0, false);
241 }
242
244}
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
UBlasMatrix< int > MatrixInt
UBlasMatrix< double > MatrixDouble
constexpr auto field_name
static MoFEMErrorCode baseFunctionsTri(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
static MoFEMErrorCode generateIndicesTriTri(const int N, int *alpha)
static MoFEMErrorCode generateIndicesEdgeTri(const int N[], int *alpha[])
const std::string fieldName