v0.16.0
Loading...
Searching...
No Matches
AnalyticalExprPython.cpp
Go to the documentation of this file.
1/**
2 * \file AnalyticalExprPython.cpp
3 * \brief Analytical expression Python bridge implementation
4 */
5
6#include <MoFEM.hpp>
7using namespace MoFEM;
8
9#include <EshelbianPlasticity.hpp>
10
11#include <regex>
12
13namespace EshelbianPlasticity {
14
15#ifdef ENABLE_PYTHON_BINDING
16
17MoFEMErrorCode
18AnalyticalExprPython::analyticalExprInit(const std::string py_file) {
19 MoFEMFunctionBegin;
20 try {
21
22 // create main module
23 auto main_module = bp::import("__main__");
24 mainNamespace = main_module.attr("__dict__");
25 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
26 // create a reference to python function
27 analyticalDispFun = mainNamespace["analytical_disp"];
28 analyticalTractionFun = mainNamespace["analytical_traction"];
29 analyticalExternalStrainFun = mainNamespace["analytical_external_strain"];
30 analyticalElasticFun = mainNamespace["analytical_elastic"];
31
32 } catch (bp::error_already_set const &) {
33 // print all other errors to stderr
34 PyErr_Print();
35 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
36 }
37 MoFEMFunctionReturn(0);
38}
39
40MoFEMErrorCode AnalyticalExprPython::evalAnalyticalDisp(
41
42 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
43 np::ndarray nx, np::ndarray ny, np::ndarray nz,
44 const std::string &block_name, np::ndarray &analytical_expr
45
46) {
47 MoFEMFunctionBegin;
48 try {
49
50 // call python function
51 analytical_expr = bp::extract<np::ndarray>(
52 analyticalDispFun(delta_t, t, x, y, z, nx, ny, nz, block_name));
53
54 } catch (bp::error_already_set const &) {
55 // print all other errors to stderr
56 PyErr_Print();
57 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
58 }
59 MoFEMFunctionReturn(0);
60}
61
62MoFEMErrorCode AnalyticalExprPython::evalAnalyticalTraction(
63
64 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
65 np::ndarray nx, np::ndarray ny, np::ndarray nz,
66 const std::string &block_name, np::ndarray &analytical_expr
67
68) {
69 MoFEMFunctionBegin;
70 try {
71
72 // call python function
73 analytical_expr = bp::extract<np::ndarray>(
74 analyticalTractionFun(delta_t, t, x, y, z, nx, ny, nz, block_name));
75
76 } catch (bp::error_already_set const &) {
77 // print all other errors to stderr
78 PyErr_Print();
79 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
80 }
81 MoFEMFunctionReturn(0);
82}
83
84MoFEMErrorCode AnalyticalExprPython::evalAnalyticalExternalStrain(
85
86 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
87 const std::string &block_name, np::ndarray &analytical_expr
88
89) {
90 MoFEMFunctionBegin;
91 try {
92
93 // call python function
94 analytical_expr = bp::extract<np::ndarray>(
95 analyticalExternalStrainFun(delta_t, t, x, y, z, block_name));
96
97 } catch (bp::error_already_set const &) {
98 // print all other errors to stderr
99 PyErr_Print();
100 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
101 }
102 MoFEMFunctionReturn(0);
103}
104
105MoFEMErrorCode AnalyticalExprPython::evalAnalyticalElastic(
106
107 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
108 const std::string &block_name, np::ndarray &analytical_expr
109
110) {
111 MoFEMFunctionBegin;
112 try {
113
114 // call python function
115 analytical_expr = bp::extract<np::ndarray>(
116 analyticalElasticFun(delta_t, t, x, y, z, block_name));
117
118 } catch (bp::error_already_set const &) {
119 // print all other errors to stderr
120 PyErr_Print();
121 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
122 }
123 MoFEMFunctionReturn(0);
124}
125
126boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
127
128namespace {
129
130np::ndarray convert_to_numpy(VectorDouble &data, int nb_gauss_pts, int id) {
131 auto dtype = np::dtype::get_builtin<double>();
132 auto size = bp::make_tuple(nb_gauss_pts);
133 auto stride = bp::make_tuple(3 * sizeof(double));
134 return (np::from_data(&data[id], dtype, size, stride, bp::object()));
135};
136
137} // namespace
138
139#endif // ENABLE_PYTHON_BINDING
140
141MatrixDouble analytical_expr_function(double delta_t, double t,
142 int nb_gauss_pts,
143 MatrixDouble &m_ref_coords,
144 MatrixDouble &m_ref_normals,
145 const std::string block_name) {
146
147#ifdef ENABLE_PYTHON_BINDING
148 if (auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
149
150 VectorDouble v_ref_coords = m_ref_coords.data();
151 VectorDouble v_ref_normals = m_ref_normals.data();
152
153 bp::list python_coords;
154 bp::list python_normals;
155
156 for (int idx = 0; idx < 3; ++idx) {
157 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
158 python_normals.append(convert_to_numpy(v_ref_normals, nb_gauss_pts, idx));
159 }
160
161 auto disp_block_name = "(.*)ANALYTICAL_DISPLACEMENT(.*)";
162 auto traction_block_name = "(.*)ANALYTICAL_TRACTION(.*)";
163
164 std::regex reg_disp_name(disp_block_name);
165 std::regex reg_traction_name(traction_block_name);
166
167 np::ndarray np_analytical_expr = np::empty(
168 bp::make_tuple(nb_gauss_pts, 3), np::dtype::get_builtin<double>());
169
170 if (std::regex_match(block_name, reg_disp_name)) {
171 CHK_MOAB_THROW(analytical_expr_ptr->evalAnalyticalDisp(
172 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
173 bp::extract<np::ndarray>(python_coords[1]),
174 bp::extract<np::ndarray>(python_coords[2]),
175 bp::extract<np::ndarray>(python_normals[0]),
176 bp::extract<np::ndarray>(python_normals[1]),
177 bp::extract<np::ndarray>(python_normals[2]),
178 block_name, np_analytical_expr),
179 "Failed analytical_disp() python call");
180 } else if (std::regex_match(block_name, reg_traction_name)) {
181 CHK_MOAB_THROW(analytical_expr_ptr->evalAnalyticalTraction(
182 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
183 bp::extract<np::ndarray>(python_coords[1]),
184 bp::extract<np::ndarray>(python_coords[2]),
185 bp::extract<np::ndarray>(python_normals[0]),
186 bp::extract<np::ndarray>(python_normals[1]),
187 bp::extract<np::ndarray>(python_normals[2]),
188 block_name, np_analytical_expr),
189 "Failed analytical_traction() python call");
190 } else {
191 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
192 "Unknown analytical block");
193 }
194
195 // check the shape of returned array
196 if (PetscUnlikely(np_analytical_expr.get_shape()[0] != nb_gauss_pts ||
197 np_analytical_expr.get_shape()[1] != 3)) {
198 CHK_THROW_MESSAGE(
199 MOFEM_DATA_INCONSISTENCY,
200 "Wrong shape of analytical expression returned from "
201 "python, expected: (" +
202 std::to_string(nb_gauss_pts) + ", 3), got: (" +
203 std::to_string(np_analytical_expr.get_shape()[0]) + ", " +
204 std::to_string(np_analytical_expr.get_shape()[1]) + ")");
205 }
206
207 double *analytical_expr_val_ptr =
208 reinterpret_cast<double *>(np_analytical_expr.get_data());
209
210 MatrixDouble v_analytical_expr;
211 MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
212 v_analytical_expr, nb_gauss_pts);
213 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
214 for (int idx = 0; idx < 3; ++idx)
215 v_analytical_expr(gg, idx) =
216 *(analytical_expr_val_ptr + (3 * gg + idx));
217 }
218 return v_analytical_expr;
219 } else {
220 // Returns empty vector
221 }
222#endif
223 return MatrixDouble();
224}
225
226VectorDouble analytical_externalstrain_function(double delta_t, double t,
227 int nb_gauss_pts,
228 MatrixDouble &m_ref_coords,
229 const std::string block_name) {
230
231#ifdef ENABLE_PYTHON_BINDING
232 if (auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
233
234 VectorDouble v_ref_coords = m_ref_coords.data();
235
236 bp::list python_coords;
237
238 for (int idx = 0; idx < 3; ++idx) {
239 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
240 }
241
242 auto externalstrain_block_name = "(.*)ANALYTICAL_EXTERNALSTRAIN(.*)";
243
244 std::regex reg_externalstrain_name(externalstrain_block_name);
245
246 np::ndarray np_analytical_expr = np::empty(
247 bp::make_tuple(nb_gauss_pts), np::dtype::get_builtin<double>());
248
249 if (std::regex_match(block_name, reg_externalstrain_name)) {
250 CHK_MOAB_THROW(analytical_expr_ptr->evalAnalyticalExternalStrain(
251 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
252 bp::extract<np::ndarray>(python_coords[1]),
253 bp::extract<np::ndarray>(python_coords[2]), block_name,
254 np_analytical_expr),
255 "Failed analytical_external_strain() python call");
256 } else {
257 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
258 "Unknown analytical block");
259 }
260
261 // check the shape of returned array
262 if (PetscUnlikely(np_analytical_expr.get_shape()[0] != nb_gauss_pts ||
263 np_analytical_expr.get_shape()[1] != 1)) {
264 CHK_THROW_MESSAGE(
265 MOFEM_DATA_INCONSISTENCY,
266 "Wrong shape of analytical expression returned from "
267 "python, expected: (" +
268 std::to_string(nb_gauss_pts) + ", 1), got: (" +
269 std::to_string(np_analytical_expr.get_shape()[0]) + ", " +
270 std::to_string(np_analytical_expr.get_shape()[1]) + ")");
271 }
272
273 double *analytical_expr_val_ptr =
274 reinterpret_cast<double *>(np_analytical_expr.get_data());
275
276 VectorDouble v_analytical_expr;
277 v_analytical_expr.resize(nb_gauss_pts, false);
278 for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
279 v_analytical_expr[gg] = *(analytical_expr_val_ptr + gg);
280
281 return v_analytical_expr;
282 } else {
283 // Returns empty vector
284 }
285#endif
286 return VectorDouble();
287}
288
289VectorDouble analytical_elastic_function(double delta_t, double t,
290 int nb_gauss_pts,
291 MatrixDouble &m_ref_coords,
292 const std::string block_name) {
293#ifdef ENABLE_PYTHON_BINDING
294 if (auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
295
296 VectorDouble v_ref_coords = m_ref_coords.data();
297
298 bp::list python_coords;
299
300 for (int idx = 0; idx < 3; ++idx) {
301 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
302 }
303
304 auto analytical_elastic_block_name = "(.*)ANALYTICAL_ELASTIC(.*)";
305
306 std::regex reg_youngs_modulus_name(analytical_elastic_block_name);
307
308 np::ndarray np_analytical_expr = np::empty(
309 bp::make_tuple(nb_gauss_pts), np::dtype::get_builtin<double>());
310
311 if (std::regex_match(block_name, reg_youngs_modulus_name)) {
312 CHK_MOAB_THROW(analytical_expr_ptr->evalAnalyticalElastic(
313 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
314 bp::extract<np::ndarray>(python_coords[1]),
315 bp::extract<np::ndarray>(python_coords[2]), block_name,
316 np_analytical_expr),
317 "Failed analytical_elastic() python call");
318 } else {
319 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
320 "Unknown analytical block");
321 }
322
323 // check the shape of returned array
324 if (PetscUnlikely(np_analytical_expr.get_shape()[0] != nb_gauss_pts ||
325 np_analytical_expr.get_shape()[1] != 1)) {
326 CHK_THROW_MESSAGE(
327 MOFEM_DATA_INCONSISTENCY,
328 "Wrong shape of analytical expression returned from "
329 "python, expected: (" +
330 std::to_string(nb_gauss_pts) + ", 1), got: (" +
331 std::to_string(np_analytical_expr.get_shape()[0]) + ", " +
332 std::to_string(np_analytical_expr.get_shape()[1]) + ")");
333 }
334
335 double *analytical_expr_val_ptr =
336 reinterpret_cast<double *>(np_analytical_expr.get_data());
337
338 VectorDouble v_analytical_expr;
339 v_analytical_expr.resize(nb_gauss_pts, false);
340 for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
341 v_analytical_expr[gg] = *(analytical_expr_val_ptr + gg);
342
343 return v_analytical_expr;
344 } else {
345 // Returns empty vector
346 }
347#endif
348 return VectorDouble();
349}
350
351} // namespace EshelbianPlasticity
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition definitions.h:592
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
EshelbianPlasticity::analytical_externalstrain_function
VectorDouble analytical_externalstrain_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_ref_coords, const std::string block_name)
EshelbianPlasticity::analytical_elastic_function
VectorDouble analytical_elastic_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_ref_coords, const std::string block_name)
EshelbianPlasticity::analytical_expr_function
MatrixDouble analytical_expr_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_ref_coords, MatrixDouble &m_ref_normals, const std::string block_name)
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition Types.hpp:68
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::GetFTensor1FromMatType
decltype(GetFTensor1FromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor1FromMatType
Definition Templates.hpp:385
analytical_expr
Definition analytical_expr.py:1
t
constexpr double t
plate stiffness
Definition plate.cpp:58
Generated by Doxygen 1.9.8 and hosted at University of Glasgow