PolyFEM
Loading...
Searching...
No Matches
q_bases.py
Go to the documentation of this file.
1# https://raw.githubusercontent.com/sympy/sympy/master/examples/advanced/fem.py
2from sympy import *
3import os
4import argparse
5from sympy.printing import ccode
6
7
8import pretty_print
9
10x, y, z = symbols('x,y,z')
11
12
13class ReferenceSimplex:
14 def __init__(self, nsd):
15 self.nsd = nsd
16 if nsd <= 3:
17 coords = symbols('x,y,z')[:nsd]
18 else:
19 coords = [Symbol("x_%d" % d) for d in range(nsd)]
20 self.coords = coords
21
22 def integrate(self, f):
23 coords = self.coords
24 nsd = self.nsd
25
26 limit = 1
27 for p in coords:
28 limit -= p
29
30 intf = f
31 for d in range(0, nsd):
32 p = coords[d]
33 limit += p
34 intf = integrate(intf, (p, 0, limit))
35 return intf
36
37
38def create_point_set(order, nsd):
39 h = Rational(1, order)
40 set = []
41
42 if nsd == 1:
43 for i in range(0, order + 1):
44 x = i * h
45 set.append([x])
46
47 if nsd == 2:
48 for i in range(0, order + 1):
49 x = i * h
50 for j in range(0, order + 1):
51 y = j * h
52 set.append((x, y))
53
54 if nsd == 3:
55 for i in range(0, order + 1):
56 x = i * h
57 for j in range(0, order + 1):
58 y = j * h
59 for k in range(0, order + 1):
60 z = k * h
61 set.append((x, y, z))
62
63 return set
64
65
66class Lagrange:
67 def __init__(self, nsd, order):
68 self.nsd = nsd
69 self.order = order
70 self.points = []
71 self.compute_basis()
72
73 def nbf(self):
74 return len(self.N)
75
76 def compute_basis(self):
77 order = self.order
78 nsd = self.nsd
79 N = []
80 self.points = create_point_set(order, nsd)
81
82 if nsd == 1:
83 Ntmpx = []
84
85 for j in range(order + 1):
86 vx = 1
87 xj = j/(order)
88 for m in range(order+1):
89 if m == j:
90 continue
91 xm = m/(order)
92 vx *= (x - xm)/(xj - xm)
93
94 Ntmpx.append(vx)
95
96 for i in range(order + 1):
97 N.append(Ntmpx[i])
98
99 elif nsd == 2:
100 Ntmpx = []
101 Ntmpy = []
102
103 for j in range(order + 1):
104 vx = 1
105 vy = 1
106 xj = 1./(order)*j
107 for m in range(order+1):
108 if m == j:
109 continue
110 xm = 1./(order)*m
111 vx *= (x - xm)/(xj - xm)
112 vy *= (y - xm)/(xj - xm)
113
114 Ntmpx.append(vx)
115 Ntmpy.append(vy)
116
117 for i in range(order + 1):
118 for j in range(order + 1):
119 N.append(Ntmpx[i]*Ntmpy[j])
120 elif nsd == 3:
121 Ntmpx = []
122 Ntmpy = []
123 Ntmpz = []
124
125 for j in range(order + 1):
126 vx = 1
127 vy = 1
128 vz = 1
129 xj = 1./(order)*j
130 for m in range(order+1):
131 if m == j:
132 continue
133 xm = 1./(order)*m
134 vx *= (x - xm)/(xj - xm)
135 vy *= (y - xm)/(xj - xm)
136 vz *= (z - xm)/(xj - xm)
137
138 Ntmpx.append(vx)
139 Ntmpy.append(vy)
140 Ntmpz.append(vz)
141
142 for i in range(order + 1):
143 for j in range(order + 1):
144 for l in range(order + 1):
145 N.append(Ntmpx[i]*Ntmpy[j]*Ntmpz[l])
146
147 self.N = N
148
149
150def parse_args():
151 parser = argparse.ArgumentParser(
152 description=__doc__,
153 formatter_class=argparse.RawDescriptionHelpFormatter)
154 parser.add_argument("output", type=str, help="path to the output folder")
155 return parser.parse_args()
156
157
158if __name__ == "__main__":
159 args = parse_args()
160 path = os.path.abspath(args.output)
161
162 dims = [1, 2, 3]
163 orders = [0, 1, 2, 3, -2]
164
165 for dim in dims:
166 namev = f"auto_q_bases_{dim}d_val"
167 namen = f"auto_q_bases_{dim}d_nodes"
168 nameg = f"auto_q_bases_{dim}d_grad"
169
170 cppv = f"#include \"{namev}.hpp\"\n\n\n"
171 cppv = cppv + \
172 "namespace polyfem {\nnamespace autogen " + "{\nnamespace " + "{\n"
173
174 cppn = f"#include \"{namen}.hpp\"\n\n\n"
175 cppn = cppn + \
176 "namespace polyfem {\nnamespace autogen " + "{\nnamespace " + "{\n"
177
178 cppg = f"#include \"{nameg}.hpp\"\n\n\n"
179 cppg = cppg + \
180 "namespace polyfem {\nnamespace autogen " + "{\nnamespace " + "{\n"
181 if dim == 3:
182 cppg = "#include <Eigen/Dense>\n#include <cassert>\n namespace polyfem {\nnamespace autogen {"
183
184 eextern = ""
185
186 hppv = "#pragma once\n\n#include <Eigen/Dense>\n#include <cassert>\n\n"
187 hppv = hppv + "namespace polyfem {\nnamespace autogen " + "{\n"
188
189 hppn = "#pragma once\n\n#include <Eigen/Dense>\n#include <cassert>\n\n"
190 hppn = hppn + "namespace polyfem {\nnamespace autogen " + "{\n"
191
192 hppg = "#pragma once\n\n#include <Eigen/Dense>\n#include <cassert>\n\n"
193 hppg = hppg + "namespace polyfem {\nnamespace autogen " + "{\n"
194
195 print(str(dim) + "D")
196 suffix = "_2d" if dim == 2 else ("_3d" if dim == 3 else "_1d")
197
198 unique_nodes = "void q_nodes" + suffix + \
199 "(const int q, Eigen::MatrixXd &val)"
200
201 unique_fun = "void q_basis_value" + suffix + \
202 "(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)"
203 dunique_fun = "void q_grad_basis_value" + suffix + \
204 "(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)"
205
206 hppn = hppn + unique_nodes + ";\n\n"
207
208 hppv = hppv + unique_fun + ";\n\n"
209 hppg = hppg + dunique_fun + ";\n\n"
210
211 unique_nodes = unique_nodes + "{\nswitch(q)" + "{\n"
212
213 unique_fun = unique_fun + "{\nswitch(q)" + "{\n"
214 dunique_fun = dunique_fun + "{\nswitch(q)" + "{\n"
215
216 if dim == 1:
217 vertices = [[0], [1]]
218 elif dim == 2:
219 vertices = [[0, 0], [1, 0], [1, 1], [0, 1]]
220 elif dim == 3:
221 vertices = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0],
222 [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]]
223
224 if dim == 1:
225 orderss = orders + [4, 5, 6]
226 else:
227 orderss = orders
228 for order in orderss:
229 print("\t-processing " + str(order))
230
231 if order == 0:
232 def fe(): return None
233 fe.nbf = lambda: 1
234
235 fe.N = [1]
236
237 if dim == 1:
238 fe.points = [[1./2]]
239 elif dim == 2:
240 fe.points = [[1./2., 1./2.]]
241 else:
242 fe.points = [[1./2., 1./2., 1./2.]]
243 elif order == -2:
244 def fe(): return None
245 if dim == 1: # no serendipity in 1d
246 fe = Lagrange(dim, 2)
247 elif dim == 2:
248 fe.points = []
249 fe.N = []
250 fe.nbf = lambda: 8
251
252 for xi_a in [-1, 1]:
253 for eta_a in [-1, 1]:
254 tmp = (1/4*(2*x*xi_a-xi_a+1))*(2*eta_a*y -
255 eta_a+1)*(2*eta_a*y+2*x*xi_a-eta_a-xi_a-1)
256 fe.N.append(tmp)
257 fe.points.append([(xi_a+1)/2, (eta_a+1)/2])
258
259 for eta_a in [-1, 1]:
260 tmp = -2*x*(x-1)*(2*eta_a*y-eta_a+1)
261 fe.N.append(tmp)
262 fe.points.append([1/2, (eta_a+1)/2])
263
264 for xi_a in [-1, 1]:
265 tmp = -2*y*(y-1)*(2*x*xi_a-xi_a+1)
266 fe.N.append(tmp)
267 fe.points.append([(xi_a+1)/2, 1/2])
268
269 assert (len(fe.points) == 8)
270 assert (len(fe.N) == 8)
271
272 elif dim == 3:
273 fe.points = []
274 fe.N = []
275 fe.nbf = lambda: 20
276
277 for xi_a in [-1, 1]:
278 for eta_a in [-1, 1]:
279 for zeta_a in [-1, 1]:
280 tmp = (1/8*(2*x*xi_a-xi_a+1))*(2*eta_a*y-eta_a+1)*(2*z*zeta_a-zeta_a+1)*(
281 2*eta_a*y+2*x*xi_a+2*z*zeta_a-eta_a-xi_a-zeta_a-2)
282 fe.N.append(tmp)
283 fe.points.append(
284 [(xi_a+1)/2, (eta_a+1)/2, (zeta_a+1)/2])
285
286 for eta_a in [-1, 1]:
287 for zeta_a in [-1, 1]:
288 tmp = -x*(x-1)*(2*eta_a*y-eta_a+1) * \
289 (2*z*zeta_a-zeta_a+1)
290 fe.N.append(tmp)
291 fe.points.append([1/2, (eta_a+1)/2, (zeta_a+1)/2])
292
293 for xi_a in [-1, 1]:
294 for zeta_a in [-1, 1]:
295 tmp = -y*(y-1)*(2*x*xi_a-xi_a+1) * \
296 (2*z*zeta_a-zeta_a+1)
297 fe.N.append(tmp)
298 fe.points.append([(xi_a+1)/2, 1/2, (zeta_a+1)/2])
299
300 for xi_a in [-1, 1]:
301 for eta_a in [-1, 1]:
302 tmp = -z*(z-1)*(2*x*xi_a-xi_a+1) * \
303 (2*eta_a*y-eta_a+1)
304 fe.N.append(tmp)
305 fe.points.append([(xi_a+1)/2, (eta_a+1)/2, 1/2])
306
307 assert (len(fe.points) == 20)
308 assert (len(fe.N) == 20)
309 else:
310 assert (False)
311 else:
312 fe = Lagrange(dim, order)
313
314 current_indices = list(range(0, len(fe.points)))
315 indices = []
316
317 if dim > 1:
318 # vertex coordinate
319 for i in range(0, 2**dim):
320 vv = vertices[i]
321
322 for ii in current_indices:
323 norm = 0
324 for dd in range(0, dim):
325 norm = norm + (vv[dd] - fe.points[ii][dd]) ** 2
326
327 if norm < 1e-10:
328 indices.append(ii)
329 current_indices.remove(ii)
330 break
331
332 # edge 0 coordinate
333 for i in range(0, abs(order) - 1):
334 for ii in current_indices:
335 if fe.points[ii][1] != 0 or (dim == 3 and fe.points[ii][2] != 0):
336 continue
337
338 if abs(fe.points[ii][0] - (i + 1) / abs(order)) < 1e-10:
339 indices.append(ii)
340 current_indices.remove(ii)
341 break
342
343 # edge 1 coordinate
344 for i in range(0, abs(order) - 1):
345 for ii in current_indices:
346 if fe.points[ii][0] != 1 or (dim == 3 and fe.points[ii][2] != 0):
347 continue
348
349 if abs(fe.points[ii][1] - (i + 1) / abs(order)) < 1e-10:
350 indices.append(ii)
351 current_indices.remove(ii)
352 break
353
354 # edge 2 coordinate
355 for i in range(0, abs(order) - 1):
356 for ii in current_indices:
357 if fe.points[ii][1] != 1 or (dim == 3 and fe.points[ii][2] != 0):
358 continue
359
360 if abs(fe.points[ii][0] - (1 - (i + 1) / abs(order))) < 1e-10:
361 indices.append(ii)
362 current_indices.remove(ii)
363 break
364
365 # edge 3 coordinate
366 for i in range(0, abs(order) - 1):
367 for ii in current_indices:
368 if fe.points[ii][0] != 0 or (dim == 3 and fe.points[ii][2] != 0):
369 continue
370
371 if abs(fe.points[ii][1] - (1 - (i + 1) / abs(order))) < 1e-10:
372 indices.append(ii)
373 current_indices.remove(ii)
374 break
375
376 if dim == 3:
377 # edge 4 coordinate
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] != 0:
381 continue
382
383 if abs(fe.points[ii][2] - (i + 1) / abs(order)) < 1e-10:
384 indices.append(ii)
385 current_indices.remove(ii)
386 break
387
388 # edge 5 coordinate
389 for i in range(0, abs(order) - 1):
390 for ii in current_indices:
391 if fe.points[ii][0] != 1 or fe.points[ii][1] != 0:
392 continue
393
394 if abs(fe.points[ii][2] - (1 - (i + 1) / abs(order))) < 1e-10:
395 indices.append(ii)
396 current_indices.remove(ii)
397 break
398
399 # edge 6 coordinate
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][1] != 1:
403 continue
404
405 if abs(fe.points[ii][2] - (1 - (i + 1) / abs(order))) < 1e-10:
406 indices.append(ii)
407 current_indices.remove(ii)
408 break
409
410 # edge 7 coordinate
411 for i in range(0, abs(order) - 1):
412 for ii in current_indices:
413 if fe.points[ii][0] != 0 or fe.points[ii][1] != 1:
414 continue
415
416 if abs(fe.points[ii][2] - (1 - (i + 1) / abs(order))) < 1e-10:
417 indices.append(ii)
418 current_indices.remove(ii)
419 break
420
421 # edge 8 coordinate
422 for i in range(0, abs(order) - 1):
423 for ii in current_indices:
424 if fe.points[ii][1] != 0 or fe.points[ii][2] != 1:
425 continue
426
427 if abs(fe.points[ii][0] - (i + 1) / abs(order)) < 1e-10:
428 indices.append(ii)
429 current_indices.remove(ii)
430 break
431
432 # edge 9 coordinate
433 for i in range(0, abs(order) - 1):
434 for ii in current_indices:
435 if fe.points[ii][0] != 1 or fe.points[ii][2] != 1:
436 continue
437
438 if abs(fe.points[ii][1] - (i + 1) / abs(order)) < 1e-10:
439 indices.append(ii)
440 current_indices.remove(ii)
441 break
442
443 # edge 10 coordinate
444 for i in range(0, abs(order) - 1):
445 for ii in current_indices:
446 if fe.points[ii][1] != 1 or fe.points[ii][2] != 1:
447 continue
448
449 if abs(fe.points[ii][0] - (1 - (i + 1) / abs(order))) < 1e-10:
450 indices.append(ii)
451 current_indices.remove(ii)
452 break
453
454 # edge 11 coordinate
455 for i in range(0, abs(order) - 1):
456 for ii in current_indices:
457 if fe.points[ii][0] != 0 or fe.points[ii][2] != 1:
458 continue
459
460 if abs(fe.points[ii][1] - (1 - (i + 1) / abs(order))) < 1e-10:
461 indices.append(ii)
462 current_indices.remove(ii)
463 break
464
465 if dim == 3:
466 nn = max(0, abs(order) - 1)
467 npts = int(nn * nn)
468
469 # side: x = 0
470 tmp = []
471 for i in range(0, npts):
472 for ii in current_indices:
473 if abs(fe.points[ii][0]) > 1e-10:
474 continue
475
476 tmp.append(ii)
477 current_indices.remove(ii)
478 break
479 tmp.sort(reverse=True)
480 indices.extend(tmp)
481
482 # side: x = 1
483 tmp = []
484 for i in range(0, npts):
485 for ii in current_indices:
486 if abs(fe.points[ii][0] - 1) > 1e-10:
487 continue
488
489 tmp.append(ii)
490 current_indices.remove(ii)
491 break
492 indices.extend(tmp)
493
494 # front: y = 0
495 for i in range(0, npts):
496 for ii in current_indices:
497 if abs(fe.points[ii][1]) > 1e-10:
498 continue
499
500 indices.append(ii)
501 current_indices.remove(ii)
502 break
503
504 # back: y = 1
505 for i in range(0, npts):
506 for ii in current_indices:
507 if abs(fe.points[ii][1]-1) > 1e-10:
508 continue
509
510 indices.append(ii)
511 current_indices.remove(ii)
512 break
513
514 # bottom: z = 0
515 for i in range(0, npts):
516 for ii in current_indices:
517 if abs(fe.points[ii][2]) > 1e-10:
518 continue
519
520 indices.append(ii)
521 current_indices.remove(ii)
522 break
523
524 # top: z = 1
525 for i in range(0, npts):
526 for ii in current_indices:
527 if abs(fe.points[ii][2]-1) > 1e-10:
528 continue
529
530 indices.append(ii)
531 current_indices.remove(ii)
532 break
533 else:
534 for i in range(0, abs(order)):
535 for ii in current_indices:
536 if abs(fe.points[ii][0] - (i / abs(order))) < 1e-10:
537 indices.append(ii)
538 current_indices.remove(ii)
539 break
540
541 # either face or volume indices, order do not matter
542 for ii in current_indices:
543 indices.append(ii)
544
545 orderN = str(order) if order >= 0 else "m"+str(-order)
546 # nodes code gen
547 nodes = "void q_" + orderN + "_nodes" + suffix + \
548 "(Eigen::MatrixXd &res) {\n res.resize(" + \
549 str(len(indices)) + ", " + str(dim) + "); res << \n"
550 unique_nodes = unique_nodes + "\tcase " + \
551 str(order) + ": " + "q_" + orderN + \
552 "_nodes" + suffix + "(val); break;\n"
553
554 eextern = eextern + \
555 f"extern \"C++\" void q_{orderN}_basis_grad_value_3d(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val);\n"
556
557 for ii in indices:
558 nodes = nodes + ccode(fe.points[ii][0]) + ("" if dim == 1 else (", " + ccode(fe.points[ii][1]) + (
559 (", " + ccode(fe.points[ii][2])) if dim == 3 else ""))) + ",\n"
560 nodes = nodes[:-2]
561 nodes = nodes + ";\n}"
562
563 # bases code gen
564 func = "void q_" + orderN + "_basis_value" + suffix + \
565 "(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &result_0)"
566 dfunc = "void q_" + orderN + "_basis_grad_value" + suffix + \
567 "(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)"
568
569 unique_fun = unique_fun + "\tcase " + \
570 str(order) + ": " + "q_" + orderN + "_basis_value" + \
571 suffix + "(local_index, uv, val); break;\n"
572 dunique_fun = dunique_fun + "\tcase " + \
573 str(order) + ": " + "q_" + orderN + "_basis_grad_value" + \
574 suffix + "(local_index, uv, val); break;\n"
575
576 # hpp = hpp + func + ";\n"
577 # hpp = hpp + dfunc + ";\n"
578
579 base = "auto x=uv.col(0).array();"
580 if dim > 1:
581 base = base + "\nauto y=uv.col(1).array();"
582 if dim == 3:
583 base = base + "\nauto z=uv.col(2).array();"
584 base = base + "\n\n"
585 dbase = base
586
587 if order == 0:
588 base = base + "result_0.resize(x.size(),1);\n"
589
590 base = base + "switch(local_index){\n"
591 dbase = dbase + \
592 "val.resize(uv.rows(), uv.cols());\n Eigen::ArrayXd result_0(uv.rows());\n" + \
593 "switch(local_index){\n"
594
595 for i in range(0, fe.nbf()):
596 real_index = indices[i]
597 # real_index = i
598
599 if dim == 3:
600 base = base + "\tcase " + str(i) + ": {" + pretty_print.C99_print(
601 simplify(fe.N[real_index])).replace(" = 1;", ".setOnes();") + "} break;\n"
602 dbase = dbase + "\tcase " + str(i) + ": {" + \
603 "{" + 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; }" \
604 "{" + pretty_print.C99_print(simplify(diff(fe.N[real_index], y))).replace(" = 0;", ".setZero();").replace(
605 " = 1;", ".setOnes();").replace(" = -1;", ".setConstant(-1);") + "val.col(1) = result_0; }"
606 dbase = dbase + "{" + pretty_print.C99_print(simplify(diff(fe.N[real_index], z))).replace(" = 0;", ".setZero();").replace(
607 " = 1;", ".setOnes();").replace(" = -1;", ".setConstant(-1);") + "val.col(2) = result_0; }"
608 elif dim == 2:
609 base = base + "\tcase " + str(i) + ": {" + pretty_print.C99_print(
610 simplify(fe.N[real_index])).replace(" = 1;", ".setOnes();") + "} break;\n"
611 dbase = dbase + "\tcase " + str(i) + ": {" + \
612 "{" + 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; }" \
613 "{" + pretty_print.C99_print(simplify(diff(fe.N[real_index], y))).replace(" = 0;", ".setZero();").replace(
614 " = 1;", ".setOnes();").replace(" = -1;", ".setConstant(-1);") + "val.col(1) = result_0; }"
615 else:
616 base = base + "\tcase " + str(i) + ": {" + pretty_print.C99_print(
617 simplify(fe.N[real_index])).replace(" = 1;", ".setOnes();") + "} break;\n"
618 dbase = dbase + "\tcase " + str(i) + ": {" + \
619 "{" + pretty_print.C99_print(simplify(diff(fe.N[real_index], x))).replace(" = 0;", ".setZero();").replace(
620 " = 1.0;", ".setOnes();").replace(" = -1.0;", ".setConstant(-1);") + "val.col(0) = result_0; }"
621
622 dbase = dbase + "} break;\n"
623
624 base = base + "\tdefault: assert(false);\n}"
625 dbase = dbase + "\tdefault: assert(false);\n}"
626
627 cppv = cppv + func + "{\n\n"
628 cppv = cppv + base + "}\n"
629
630 cppg = cppg + dfunc + "{\n\n"
631 cppg = cppg + dbase + "}\n\n"
632 cppn = cppn + nodes + "\n\n"
633
634 if dim == 3:
635 with open(os.path.join(path, f"{nameg}_{order}.cpp"), "w") as file:
636 file.write(cppg+"}}")
637 cppg = "#include <Eigen/Dense>\n#include <cassert>\n namespace polyfem {\nnamespace autogen {"
638
639 if dim == 3:
640 cppg = ""
641 unique_nodes = unique_nodes + "\tdefault: assert(false);\n}}"
642
643 unique_fun = unique_fun + "\tdefault: assert(false);\n}}"
644 dunique_fun = dunique_fun + "\tdefault: assert(false);\n}}"
645
646 cppv = cppv + "}\n\n"
647 cppn = cppn + "}\n\n"
648 if dim != 3:
649 cppg = cppg + "}\n\n"
650
651 cppn = cppn + unique_nodes + "\n}}\n"
652 cppv = cppv + unique_fun + "\n}}\n"
653 cppg = cppg + dunique_fun + "\n}}\n"
654 hppv = hppv + "\n}}\n"
655 hppn = hppn + "\n}}\n"
656 hppg = hppg + "\n}}\n"
657
658 if dim == 3:
659 tcppg = f"#include \"{nameg}.hpp\"\n\n\n"
660 tcppg = tcppg + "namespace polyfem {\nnamespace autogen {\n"
661 tcppg = tcppg + eextern + "\n"
662 cppg = tcppg+cppg
663
664 print("saving...")
665 with open(os.path.join(path, f"{namev}.cpp"), "w") as file:
666 file.write(cppv)
667 with open(os.path.join(path, f"{namen}.cpp"), "w") as file:
668 file.write(cppn)
669 with open(os.path.join(path, f"{nameg}.cpp"), "w") as file:
670 file.write(cppg)
671
672 with open(os.path.join(path, f"{namev}.hpp"), "w") as file:
673 file.write(hppv)
674 with open(os.path.join(path, f"{namen}.hpp"), "w") as file:
675 file.write(hppn)
676 with open(os.path.join(path, f"{nameg}.hpp"), "w") as file:
677 file.write(hppg)
678
679 hpp = "#pragma once\n\n#include <Eigen/Dense>\n#include <cassert>\n\n"
680 for dim in dims:
681 hpp = hpp + f"#include \"auto_q_bases_{dim}d_val.hpp\"\n"
682 hpp = hpp + f"#include \"auto_q_bases_{dim}d_nodes.hpp\"\n"
683 hpp = hpp + f"#include \"auto_q_bases_{dim}d_grad.hpp\"\n"
684 hpp = hpp + "\n\nnamespace polyfem {\nnamespace autogen " + "{\n"
685 hpp = hpp + "\nstatic const int MAX_Q_BASES = " + str(max(orders)) + ";\n"
686 hpp = hpp + "\n}}\n"
687
688 print("saving...")
689 with open(os.path.join(path, "auto_q_bases.hpp"), "w") as file:
690 file.write(hpp)
691
692 print("done!")
q_bases.Lagrange
Definition q_bases.py:66
q_bases.Lagrange.__init__
__init__(self, nsd, order)
Definition q_bases.py:67
q_bases.Lagrange.compute_basis
compute_basis(self)
Definition q_bases.py:76
q_bases.Lagrange.N
N
Definition q_bases.py:74
q_bases.Lagrange.order
order
Definition q_bases.py:69
q_bases.Lagrange.points
points
Definition q_bases.py:70
q_bases.Lagrange.nsd
nsd
Definition q_bases.py:68
q_bases.ReferenceSimplex
Definition q_bases.py:13
q_bases.ReferenceSimplex.__init__
__init__(self, nsd)
Definition q_bases.py:14
q_bases.ReferenceSimplex.integrate
integrate(self, f)
Definition q_bases.py:22
q_bases.ReferenceSimplex.nsd
nsd
Definition q_bases.py:15
q_bases.ReferenceSimplex.coords
coords
Definition q_bases.py:20
pretty_print.C99_print
C99_print(expr)
Definition pretty_print.py:8
q_bases.create_point_set
create_point_set(order, nsd)
Definition q_bases.py:38
q_bases.parse_args
parse_args()
Definition q_bases.py:150
q_bases.fe
fe
Definition q_bases.py:246
q_bases.nbf
nbf
Definition q_bases.py:233