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
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
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
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\nnamespace polyfem {\nnamespace autogen {\nnamespace {\n"
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 # Generate two functions:
565 # - "func" to eval basis value
566 # - "dfunc" to eval basis gradient.
567 # Both function evaluates quadrature points in batch by dispatching the "single" basis
568 # kernel inside a for loop. In this script the single kernel related codegen variable
569 # is denoted with scalar_ prefix.
570 func = "void q_" + orderN + "_basis_value" + suffix + \
571 "(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &result_0)"
572 dfunc = "void q_" + orderN + "_basis_grad_value" + suffix + \
573 "(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)"
574 scalar_func_name = f"q_{orderN}_basis_value{suffix}_single"
575 scalar_dfunc_name = f"q_{orderN}_basis_grad_value{suffix}_single"
576
577 unique_fun = unique_fun + "\tcase " + \
578 str(order) + ": " + "q_" + orderN + "_basis_value" + \
579 suffix + "(local_index, uv, val); break;\n"
580 dunique_fun = dunique_fun + "\tcase " + \
581 str(order) + ": " + "q_" + orderN + "_basis_grad_value" + \
582 suffix + "(local_index, uv, val); break;\n"
583
584 # hpp = hpp + func + ";\n"
585 # hpp = hpp + dfunc + ";\n"
586
587 # Single basis kernel.
588 base = ""
589 dbase = ""
590 # Batch basis kernel. Basically a switch dispatch + for loop.
591 base_cases = "switch(local_index){\n"
592 dbase_cases = "switch(local_index){\n"
593
594 for i in range(0, fe.nbf()):
595 real_index = indices[i]
596 value_name = f"{scalar_func_name}_{i}"
597 grad_name = f"{scalar_dfunc_name}_{i}"
598 basis = simplify(fe.N[real_index])
599
600 base = base + pretty_print.C99_print_scalar_value_function(value_name, basis, dim)
601 dbase = dbase + pretty_print.C99_print_scalar_gradient_function(grad_name, basis, dim)
602 base_cases = base_cases + pretty_print.C99_print_scalar_value_case(i, value_name, dim)
603 dbase_cases = dbase_cases + pretty_print.C99_print_scalar_gradient_case(i, grad_name, dim)
604
605 base_cases = base_cases + "\tdefault: assert(false);\n}"
606 dbase_cases = dbase_cases + "\tdefault: assert(false);\n}"
607
608 cppv = cppv + base + "\n\n"
609 cppv = cppv + func + "{\n"
610 cppv = cppv + "result_0.resize(uv.rows(), 1);\n"
611 cppv = cppv + base_cases + "\n}\n"
612
613 cppg = cppg + dbase + "\n\n"
614 if dim == 3:
615 cppg = cppg + "}\n\n"
616 cppg = cppg + dfunc + "{\n"
617 cppg = cppg + f"val.resize(uv.rows(), {dim});\n"
618 cppg = cppg + f"double gradient[{dim}];\n"
619 cppg = cppg + dbase_cases + "\n}\n\n"
620 cppn = cppn + nodes + "\n\n"
621
622 if dim == 3:
623 with open(os.path.join(path, f"{nameg}_{order}.cpp"), "w") as file:
624 file.write(cppg+"}}")
625 cppg = "#include <Eigen/Dense>\n#include <cassert>\n\nnamespace polyfem {\nnamespace autogen {\nnamespace {\n"
626
627 if dim == 3:
628 cppg = ""
629 unique_nodes = unique_nodes + "\tdefault: assert(false);\n}}"
630
631 unique_fun = unique_fun + "\tdefault: assert(false);\n}}"
632 dunique_fun = dunique_fun + "\tdefault: assert(false);\n}}"
633
634 cppv = cppv + "}\n\n"
635 cppn = cppn + "}\n\n"
636 if dim != 3:
637 cppg = cppg + "}\n\n"
638
639 cppn = cppn + unique_nodes + "\n}}\n"
640 cppv = cppv + unique_fun + "\n}}\n"
641 cppg = cppg + dunique_fun + "\n}}\n"
642 hppv = hppv + "\n}}\n"
643 hppn = hppn + "\n}}\n"
644 hppg = hppg + "\n}}\n"
645
646 if dim == 3:
647 tcppg = f"#include \"{nameg}.hpp\"\n\n\n"
648 tcppg = tcppg + "namespace polyfem {\nnamespace autogen {\n"
649 tcppg = tcppg + eextern + "\n"
650 cppg = tcppg+cppg
651
652 print("saving...")
653 with open(os.path.join(path, f"{namev}.cpp"), "w") as file:
654 file.write(cppv)
655 with open(os.path.join(path, f"{namen}.cpp"), "w") as file:
656 file.write(cppn)
657 with open(os.path.join(path, f"{nameg}.cpp"), "w") as file:
658 file.write(cppg)
659
660 with open(os.path.join(path, f"{namev}.hpp"), "w") as file:
661 file.write(hppv)
662 with open(os.path.join(path, f"{namen}.hpp"), "w") as file:
663 file.write(hppn)
664 with open(os.path.join(path, f"{nameg}.hpp"), "w") as file:
665 file.write(hppg)
666
667 hpp = "#pragma once\n\n#include <Eigen/Dense>\n#include <cassert>\n\n"
668 for dim in dims:
669 hpp = hpp + f"#include \"auto_q_bases_{dim}d_val.hpp\"\n"
670 hpp = hpp + f"#include \"auto_q_bases_{dim}d_nodes.hpp\"\n"
671 hpp = hpp + f"#include \"auto_q_bases_{dim}d_grad.hpp\"\n"
672 hpp = hpp + "\n\nnamespace polyfem {\nnamespace autogen " + "{\n"
673 hpp = hpp + "\nstatic const int MAX_Q_BASES = " + str(max(orders)) + ";\n"
674 hpp = hpp + "\n}}\n"
675
676 print("saving...")
677 with open(os.path.join(path, "auto_q_bases.hpp"), "w") as file:
678 file.write(hpp)
679
680 print("done!")
__init__(self, nsd, order)
Definition q_bases.py:67
compute_basis(self)
Definition q_bases.py:76
C99_print_scalar_gradient_function(function_name, expr, dim)
C99_print_scalar_value_function(function_name, expr, dim)
C99_print_scalar_value_case(local_index, function_name, dim)
C99_print_scalar_gradient_case(local_index, function_name, dim)
create_point_set(order, nsd)
Definition q_bases.py:38
parse_args()
Definition q_bases.py:150