5from sympy.printing
import ccode
10x, y, z = symbols(
'x,y,z')
17 coords = symbols(
'x,y,z')[:nsd]
19 coords = [Symbol(
"x_%d" % d)
for d
in range(nsd)]
31 for d
in range(0, nsd):
39 h = Rational(1, order)
43 for i
in range(0, order + 1):
45 for j
in range(0, order + 1):
50 for i
in range(0, order + 1):
52 for j
in range(0, order + 1):
54 for k
in range(0, order + 1):
81 for j
in range(order + 1):
85 for m
in range(order+1):
89 vx *= (x - xm)/(xj - xm)
90 vy *= (y - xm)/(xj - xm)
95 for i
in range(order + 1):
96 for j
in range(order + 1):
97 N.append(Ntmpx[i]*Ntmpy[j])
103 for j
in range(order + 1):
108 for m
in range(order+1):
112 vx *= (x - xm)/(xj - xm)
113 vy *= (y - xm)/(xj - xm)
114 vz *= (z - xm)/(xj - xm)
120 for i
in range(order + 1):
121 for j
in range(order + 1):
122 for l
in range(order + 1):
123 N.append(Ntmpx[i]*Ntmpy[j]*Ntmpz[l])
129 parser = argparse.ArgumentParser(
131 formatter_class=argparse.RawDescriptionHelpFormatter)
132 parser.add_argument(
"output", type=str, help=
"path to the output folder")
133 return parser.parse_args()
136if __name__ ==
"__main__":
138 path = os.path.abspath(args.output)
141 orders = [0, 1, 2, 3, -2]
144 namev = f
"auto_q_bases_{dim}d_val"
145 namen = f
"auto_q_bases_{dim}d_nodes"
146 nameg = f
"auto_q_bases_{dim}d_grad"
148 cppv = f
"#include \"{namev}.hpp\"\n\n\n"
150 "namespace polyfem {\nnamespace autogen " +
"{\nnamespace " +
"{\n"
152 cppn = f
"#include \"{namen}.hpp\"\n\n\n"
154 "namespace polyfem {\nnamespace autogen " +
"{\nnamespace " +
"{\n"
156 cppg = f
"#include \"{nameg}.hpp\"\n\n\n"
158 "namespace polyfem {\nnamespace autogen " +
"{\nnamespace " +
"{\n"
160 cppg =
"#include <Eigen/Dense>\n#include <cassert>\n namespace polyfem {\nnamespace autogen {"
164 hppv =
"#pragma once\n\n#include <Eigen/Dense>\n#include <cassert>\n\n"
165 hppv = hppv +
"namespace polyfem {\nnamespace autogen " +
"{\n"
167 hppn =
"#pragma once\n\n#include <Eigen/Dense>\n#include <cassert>\n\n"
168 hppn = hppn +
"namespace polyfem {\nnamespace autogen " +
"{\n"
170 hppg =
"#pragma once\n\n#include <Eigen/Dense>\n#include <cassert>\n\n"
171 hppg = hppg +
"namespace polyfem {\nnamespace autogen " +
"{\n"
173 print(str(dim) +
"D")
174 suffix =
"_2d" if dim == 2
else "_3d"
176 unique_nodes =
"void q_nodes" + suffix + \
177 "(const int q, Eigen::MatrixXd &val)"
179 unique_fun =
"void q_basis_value" + suffix + \
180 "(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)"
181 dunique_fun =
"void q_grad_basis_value" + suffix + \
182 "(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)"
184 hppn = hppn + unique_nodes +
";\n\n"
186 hppv = hppv + unique_fun +
";\n\n"
187 hppg = hppg + dunique_fun +
";\n\n"
189 unique_nodes = unique_nodes +
"{\nswitch(q)" +
"{\n"
191 unique_fun = unique_fun +
"{\nswitch(q)" +
"{\n"
192 dunique_fun = dunique_fun +
"{\nswitch(q)" +
"{\n"
195 vertices = [[0, 0], [1, 0], [1, 1], [0, 1]]
197 vertices = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0],
198 [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]]
201 print(
"\t-processing " + str(order))
204 def fe():
return None
210 fe.points = [[1./2., 1./2.]]
212 fe.points = [[1./2., 1./2., 1./2.]]
214 def fe():
return None
221 for eta_a
in [-1, 1]:
222 tmp = (1/4*(2*x*xi_a-xi_a+1))*(2*eta_a*y -
223 eta_a+1)*(2*eta_a*y+2*x*xi_a-eta_a-xi_a-1)
225 fe.points.append([(xi_a+1)/2, (eta_a+1)/2])
227 for eta_a
in [-1, 1]:
228 tmp = -2*x*(x-1)*(2*eta_a*y-eta_a+1)
230 fe.points.append([1/2, (eta_a+1)/2])
233 tmp = -2*y*(y-1)*(2*x*xi_a-xi_a+1)
235 fe.points.append([(xi_a+1)/2, 1/2])
237 assert (len(fe.points) == 8)
238 assert (len(fe.N) == 8)
246 for eta_a
in [-1, 1]:
247 for zeta_a
in [-1, 1]:
248 tmp = (1/8*(2*x*xi_a-xi_a+1))*(2*eta_a*y-eta_a+1)*(2*z*zeta_a-zeta_a+1)*(
249 2*eta_a*y+2*x*xi_a+2*z*zeta_a-eta_a-xi_a-zeta_a-2)
252 [(xi_a+1)/2, (eta_a+1)/2, (zeta_a+1)/2])
254 for eta_a
in [-1, 1]:
255 for zeta_a
in [-1, 1]:
256 tmp = -x*(x-1)*(2*eta_a*y-eta_a+1) * \
257 (2*z*zeta_a-zeta_a+1)
259 fe.points.append([1/2, (eta_a+1)/2, (zeta_a+1)/2])
262 for zeta_a
in [-1, 1]:
263 tmp = -y*(y-1)*(2*x*xi_a-xi_a+1) * \
264 (2*z*zeta_a-zeta_a+1)
266 fe.points.append([(xi_a+1)/2, 1/2, (zeta_a+1)/2])
269 for eta_a
in [-1, 1]:
270 tmp = -z*(z-1)*(2*x*xi_a-xi_a+1) * \
273 fe.points.append([(xi_a+1)/2, (eta_a+1)/2, 1/2])
275 assert (len(fe.points) == 20)
276 assert (len(fe.N) == 20)
283 current_indices = list(range(0, len(fe.points)))
287 for i
in range(0, 4*(dim-1)):
289 for ii
in current_indices:
291 for dd
in range(0, dim):
292 norm = norm + (vv[dd] - fe.points[ii][dd]) ** 2
296 current_indices.remove(ii)
300 for i
in range(0, abs(order) - 1):
301 for ii
in current_indices:
302 if fe.points[ii][1] != 0
or (dim == 3
and fe.points[ii][2] != 0):
305 if abs(fe.points[ii][0] - (i + 1) / abs(order)) < 1e-10:
307 current_indices.remove(ii)
311 for i
in range(0, abs(order) - 1):
312 for ii
in current_indices:
313 if fe.points[ii][0] != 1
or (dim == 3
and fe.points[ii][2] != 0):
316 if abs(fe.points[ii][1] - (i + 1) / abs(order)) < 1e-10:
318 current_indices.remove(ii)
322 for i
in range(0, abs(order) - 1):
323 for ii
in current_indices:
324 if fe.points[ii][1] != 1
or (dim == 3
and fe.points[ii][2] != 0):
327 if abs(fe.points[ii][0] - (1 - (i + 1) / abs(order))) < 1e-10:
329 current_indices.remove(ii)
333 for i
in range(0, abs(order) - 1):
334 for ii
in current_indices:
335 if fe.points[ii][0] != 0
or (dim == 3
and fe.points[ii][2] != 0):
338 if abs(fe.points[ii][1] - (1 - (i + 1) / abs(order))) < 1e-10:
340 current_indices.remove(ii)
345 for i
in range(0, abs(order) - 1):
346 for ii
in current_indices:
347 if fe.points[ii][0] != 0
or fe.points[ii][1] != 0:
350 if abs(fe.points[ii][2] - (i + 1) / abs(order)) < 1e-10:
352 current_indices.remove(ii)
356 for i
in range(0, abs(order) - 1):
357 for ii
in current_indices:
358 if fe.points[ii][0] != 1
or fe.points[ii][1] != 0:
361 if abs(fe.points[ii][2] - (1 - (i + 1) / abs(order))) < 1e-10:
363 current_indices.remove(ii)
367 for i
in range(0, abs(order) - 1):
368 for ii
in current_indices:
369 if fe.points[ii][0] != 1
or fe.points[ii][1] != 1:
372 if abs(fe.points[ii][2] - (1 - (i + 1) / abs(order))) < 1e-10:
374 current_indices.remove(ii)
378 for i
in range(0, abs(order) - 1):
379 for ii
in current_indices:
380 if fe.points[ii][0] != 0
or fe.points[ii][1] != 1:
383 if abs(fe.points[ii][2] - (1 - (i + 1) / abs(order))) < 1e-10:
385 current_indices.remove(ii)
389 for i
in range(0, abs(order) - 1):
390 for ii
in current_indices:
391 if fe.points[ii][1] != 0
or fe.points[ii][2] != 1:
394 if abs(fe.points[ii][0] - (i + 1) / abs(order)) < 1e-10:
396 current_indices.remove(ii)
400 for i
in range(0, abs(order) - 1):
401 for ii
in current_indices:
402 if fe.points[ii][0] != 1
or fe.points[ii][2] != 1:
405 if abs(fe.points[ii][1] - (i + 1) / abs(order)) < 1e-10:
407 current_indices.remove(ii)
411 for i
in range(0, abs(order) - 1):
412 for ii
in current_indices:
413 if fe.points[ii][1] != 1
or fe.points[ii][2] != 1:
416 if abs(fe.points[ii][0] - (1 - (i + 1) / abs(order))) < 1e-10:
418 current_indices.remove(ii)
422 for i
in range(0, abs(order) - 1):
423 for ii
in current_indices:
424 if fe.points[ii][0] != 0
or fe.points[ii][2] != 1:
427 if abs(fe.points[ii][1] - (1 - (i + 1) / abs(order))) < 1e-10:
429 current_indices.remove(ii)
433 nn = max(0, abs(order) - 1)
438 for i
in range(0, npts):
439 for ii
in current_indices:
440 if abs(fe.points[ii][0]) > 1e-10:
444 current_indices.remove(ii)
446 tmp.sort(reverse=
True)
451 for i
in range(0, npts):
452 for ii
in current_indices:
453 if abs(fe.points[ii][0] - 1) > 1e-10:
457 current_indices.remove(ii)
462 for i
in range(0, npts):
463 for ii
in current_indices:
464 if abs(fe.points[ii][1]) > 1e-10:
468 current_indices.remove(ii)
472 for i
in range(0, npts):
473 for ii
in current_indices:
474 if abs(fe.points[ii][1]-1) > 1e-10:
478 current_indices.remove(ii)
482 for i
in range(0, npts):
483 for ii
in current_indices:
484 if abs(fe.points[ii][2]) > 1e-10:
488 current_indices.remove(ii)
492 for i
in range(0, npts):
493 for ii
in current_indices:
494 if abs(fe.points[ii][2]-1) > 1e-10:
498 current_indices.remove(ii)
502 for ii
in current_indices:
505 orderN = str(order)
if order >= 0
else "m"+str(-order)
507 nodes =
"void q_" + orderN +
"_nodes" + suffix + \
508 "(Eigen::MatrixXd &res) {\n res.resize(" + \
509 str(len(indices)) +
", " + str(dim) +
"); res << \n"
510 unique_nodes = unique_nodes +
"\tcase " + \
511 str(order) +
": " +
"q_" + orderN + \
512 "_nodes" + suffix +
"(val); break;\n"
514 eextern = eextern + \
515 f
"extern \"C++\" void q_{orderN}_basis_grad_value_3d(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val);\n"
518 nodes = nodes + ccode(fe.points[ii][0]) +
", " + ccode(fe.points[ii][1]) + (
519 (
", " + ccode(fe.points[ii][2]))
if dim == 3
else "") +
",\n"
521 nodes = nodes +
";\n}"
524 func =
"void q_" + orderN +
"_basis_value" + suffix + \
525 "(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &result_0)"
526 dfunc =
"void q_" + orderN +
"_basis_grad_value" + suffix + \
527 "(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)"
529 unique_fun = unique_fun +
"\tcase " + \
530 str(order) +
": " +
"q_" + orderN +
"_basis_value" + \
531 suffix +
"(local_index, uv, val); break;\n"
532 dunique_fun = dunique_fun +
"\tcase " + \
533 str(order) +
": " +
"q_" + orderN +
"_basis_grad_value" + \
534 suffix +
"(local_index, uv, val); break;\n"
539 base =
"auto x=uv.col(0).array();\nauto y=uv.col(1).array();"
541 base = base +
"\nauto z=uv.col(2).array();"
546 base = base +
"result_0.resize(x.size(),1);\n"
548 base = base +
"switch(local_index){\n"
550 "val.resize(uv.rows(), uv.cols());\n Eigen::ArrayXd result_0(uv.rows());\n" + \
551 "switch(local_index){\n"
553 for i
in range(0, fe.nbf()):
554 real_index = indices[i]
559 simplify(fe.N[real_index])).replace(
" = 1;",
".setOnes();") +
"} break;\n"
560 dbase = dbase +
"\tcase " + str(i) +
": {" + \
561 "{" +
pretty_print.C99_print(simplify(diff(fe.N[real_index], x))).replace(
" = 0;",
".setZero();").replace(
" = 1;",
".setOnes();").replace(
" = -1;",
".setConstant(-1);") +
"val.col(0) = result_0; }" \
562 "{" +
pretty_print.C99_print(simplify(diff(fe.N[real_index], y))).replace(
" = 0;",
".setZero();").replace(
563 " = 1;",
".setOnes();").replace(
" = -1;",
".setConstant(-1);") +
"val.col(1) = result_0; }"
564 dbase = dbase +
"{" +
pretty_print.C99_print(simplify(diff(fe.N[real_index], z))).replace(
" = 0;",
".setZero();").replace(
565 " = 1;",
".setOnes();").replace(
" = -1;",
".setConstant(-1);") +
"val.col(2) = result_0; }"
568 simplify(fe.N[real_index])).replace(
" = 1;",
".setOnes();") +
"} break;\n"
569 dbase = dbase +
"\tcase " + str(i) +
": {" + \
570 "{" +
pretty_print.C99_print(simplify(diff(fe.N[real_index], x))).replace(
" = 0;",
".setZero();").replace(
" = 1;",
".setOnes();").replace(
" = -1;",
".setConstant(-1);") +
"val.col(0) = result_0; }" \
571 "{" +
pretty_print.C99_print(simplify(diff(fe.N[real_index], y))).replace(
" = 0;",
".setZero();").replace(
572 " = 1;",
".setOnes();").replace(
" = -1;",
".setConstant(-1);") +
"val.col(1) = result_0; }"
574 dbase = dbase +
"} break;\n"
576 base = base +
"\tdefault: assert(false);\n}"
577 dbase = dbase +
"\tdefault: assert(false);\n}"
579 cppv = cppv + func +
"{\n\n"
580 cppv = cppv + base +
"}\n"
582 cppg = cppg + dfunc +
"{\n\n"
583 cppg = cppg + dbase +
"}\n\n"
584 cppn = cppn + nodes +
"\n\n"
587 with open(os.path.join(path, f
"{nameg}_{order}.cpp"),
"w")
as file:
588 file.write(cppg+
"}}")
589 cppg =
"#include <Eigen/Dense>\n#include <cassert>\n namespace polyfem {\nnamespace autogen {"
593 unique_nodes = unique_nodes +
"\tdefault: assert(false);\n}}"
595 unique_fun = unique_fun +
"\tdefault: assert(false);\n}}"
596 dunique_fun = dunique_fun +
"\tdefault: assert(false);\n}}"
598 cppv = cppv +
"}\n\n"
599 cppn = cppn +
"}\n\n"
601 cppg = cppg +
"}\n\n"
603 cppn = cppn + unique_nodes +
"\n}}\n"
604 cppv = cppv + unique_fun +
"\n}}\n"
605 cppg = cppg + dunique_fun +
"\n}}\n"
606 hppv = hppv +
"\n}}\n"
607 hppn = hppn +
"\n}}\n"
608 hppg = hppg +
"\n}}\n"
611 tcppg = f
"#include \"{nameg}.hpp\"\n\n\n"
612 tcppg = tcppg +
"namespace polyfem {\nnamespace autogen {\n"
613 tcppg = tcppg + eextern +
"\n"
617 with open(os.path.join(path, f
"{namev}.cpp"),
"w")
as file:
619 with open(os.path.join(path, f
"{namen}.cpp"),
"w")
as file:
621 with open(os.path.join(path, f
"{nameg}.cpp"),
"w")
as file:
624 with open(os.path.join(path, f
"{namev}.hpp"),
"w")
as file:
626 with open(os.path.join(path, f
"{namen}.hpp"),
"w")
as file:
628 with open(os.path.join(path, f
"{nameg}.hpp"),
"w")
as file:
631 hpp =
"#pragma once\n\n#include <Eigen/Dense>\n#include <cassert>\n\n"
633 hpp = hpp + f
"#include \"auto_q_bases_{dim}d_val.hpp\"\n"
634 hpp = hpp + f
"#include \"auto_q_bases_{dim}d_nodes.hpp\"\n"
635 hpp = hpp + f
"#include \"auto_q_bases_{dim}d_grad.hpp\"\n"
636 hpp = hpp +
"\n\nnamespace polyfem {\nnamespace autogen " +
"{\n"
637 hpp = hpp +
"\nstatic const int MAX_Q_BASES = " + str(max(orders)) +
";\n"
641 with open(os.path.join(path,
"auto_q_bases.hpp"),
"w")
as file:
__init__(self, nsd, order)
create_point_set(order, nsd)