Peano 4
Loading...
Searching...
No Matches
LagrangeBasisWithDiagonalMassMatrix.py
Go to the documentation of this file.
1# This file is part of the ExaHyPE2 project. For conditions of distribution and
2# use, please see the copyright notice at www.peano-framework.org
3
4
5#import mpmath as mp
6#import numpy as np
7#from sympy.integrals.quadrature import gauss_legendre, gauss_lobatto # quadrature_points by default in [-1,1]
8
9
10import scipy
11from scipy.integrate import quad
12
13from abc import abstractmethod
14
15from .LagrangeBasis import LagrangeBasis
16from .LagrangeBasis import render_tensor_1
17from .LagrangeBasis import render_tensor_2
18from .LagrangeBasis import render_tensor_3
19
20from numpy import number
21
22
23
24
26 """
27
28 Lagrange Basis with a diagonal mass matrix
29
30
31 Class representing an abstract Lagrange basis which is used as a factory
32 mechanism to produce vectors and matrices that we need within DG solvers.
33 Never use this class directly, but always a subclass. As with all basis
34 classes, the most relevant
35 routine for a user is init_dictionary_with_default_parameters() which is
36 a routine to befill a map with replacement rules, i.e. with the precomputed
37 matrices and vectors. The class assumes that the mass matrix is a diagonal
38 matrix.
39 This is not true in general, i.e. holds only for Gauss-Legendre sampling
40 points, but you might want to work with lumped matrices, and then this
41 class is the right starting point.
42
43 Further to the initialisation, some using classes use quadrature_points_and_weights().
44 The exahype2.solvers.rkdg.RungeKuttaDG solver for example uses it to
45 initialise the mapping of the plotter properly.
46
47
48 ## Generated vectors and matrices
49
50 I decided to "outsource" the documentation of the individual fields to
51 class attributes of their own. This way, they are visible both in the
52 doxygen documentation and I can pipe them into the generated C++ code.
53
54 - QuadraturePoints1d
55 - MassMatrixDiagonal1d
56 - StiffnessOperator
57 - RestrictionMatrix1d
58 - InterpolationMatrix1d
59
60
61 ## Realisation
62
63 The first version of this utility class has been written by Dominic E.
64 Charrier. It relied on NumPy and SymPy in several places. The current
65 version is stripped down and solely uses SciPy.
66
67
68 """
69
70 def __init__(self,polynomial_order):
71 super(LagrangeBasisWithDiagonalMassMatrix,self).__init__(polynomial_order)
72
73
74 def __str__(self):
75 d = {}
77 return "<{}.{} object>: {}".format(self.__class__.__module__,self.__class__.__name__,d)
78
79
80 __repr__ = __str__
81
82
83 QuadraturePoints1dDocumentation = """
84 The quadrature points are an array of size [p+1] if p is the order of the
85 chosen polyomial. The double array defines the spacing of the quadrature points
86 over the unit interval. See the helper routine exahype2::dg::getQuadraturePoint()
87 for an example how the points are used to identify the actual location in
88 space of the quadrature points.
89
90 As we work in a tensor product world, the 2d and 3d quadrature points can be
91 simply constructed by multiplying the respective coordinate system indices.
92"""
93
94
95 QuadratureWeights1dDocumentation = """
96 The quadrature weights are an array of size [p+1] if p is the order of the
97 chosen polyomial. The double array defines the integral over the polynomial
98 over the unit interval. We use Lagrangian bases here, i.e. only one polynomial
99 of the basis is 1 in a quadrature point, while all the others are 0. However,
100 the fact that it is 1 in the point, does not tell us anything about the
101 integral over this particular polynomial. That's the weight associated with
102 it.
103
104 As we work in a tensor product world, the 2d and 3d weights can be
105 simply constructed by multiplying the respective coordinate system indices.
106"""
107
108
109 MassMatrix1dDocumentation = """
110 Computes the (reference) element mass matrix' diagonal. The mass matrix is defined
111 via
112
113 @f$ \int _{(0,1)}(\phi _i, \phi _j) dx, i,j = 1 ... N > 0 @f$
114
115 where @f$ \phi _i,\phi _j @f$ are Lagrange basis functions associated with
116 the quadrature within the unit-interval. They are polynomials of maximum order N-1.
117
118 Our implementation is quite some overkill, as I assemble the whole mass matrix
119 and then extract only the diagonal. If we use Legendre nodes as support points,
120 this diagonal is an exact evaluation, i.e. the mass matrix is diagonal. If
121 we use Lobatto nodes, it is an approximation.
122
123 Please note that this is a 1d mass matrix: If we have a polynomial of order p,
124 we have @f$ (p+1)^d @f$ dofs in total per cell. My matrix has only the dimensions
125 @f$ (p+1) \times (p+1) @f$ so it effectively only models how the dofs along one
126 horizontal or vertical line within the cell do interact. This is totally sufficient
127 if we know that we have a tensor product approach of the shape functions and a
128 diagonal mass matrix.
129
130 @todo I don't know if we need renormalisation for Lobatto nodes, i.e. if we
131 should ensure that the row sum is preserved
132
133 This is an operation which maps nodal data into a weak space, i.e. the outcome has
134 to be scaled with h to be correct. It stems from a 1d integral. If you want to
135 compute the real mass matrix output, you have to multiply the outcome along each
136 coordinate axis.
137"""
138
139
140 StiffnessOperator1dDocumentation = """
141 The stiffness operator is a [p+1]x[p+1] matrix which is constructed in
142 __compute_stiffness_operator().
143
144 The matrix holds
145
146 @f$ K_{ij} = <\partial _x \phi _i, \phi _j> = w_j (\partial _x \phi _j) (x_i) @f$
147
148 where @f$ \phi _i, \phi _j @f$ are Lagrange basis functions associated with
149 the support (quadrature) points @f$ x_i @f$ and @f$ x_j @f$ within the unit
150 interval.
151
152 If you have a vector of p+1 entries, it defines a polynomial (aka function f) over the unit
153 interval. If you multiply this vector if @f$ K @f$, you get the outcome of
154
155 @f$ \int _{[0.1]} \Big( \partial _x f, \phi _j \Big) dx @f$
156
157 for the different tests @f$ \phi _j @f$. If your interval is not the unit
158 one but has side length $h$, you have to multiply the outcome with $h$ as
159 you integrate. At the same time, the derivative however has to be rescaled
160 with @f$ h^{-1} @f$ and the two h-recalibrations thus cancel out.
161 If you have a 2d or 3d function space, you still have to multiple with these
162 shape functions (see discussion around mass matrix), and you will get an
163 additional @f$ h^{d-1} @f$ scaling of the outcome.
164
165 We note that the outcome is exact for both Lobatto and Legendre nodes due to
166 the reduced degree in one scalar product operand.
167"""
168
169
170 BasisFunctionValuesLeft1dDocumentation = """
171 This is a vector of p+1 entries which tells you for a 1d interval
172 what the function on the left side is. Assume all Lagrange polynomials
173 are defined over the unit interval. The first entry in this matrix
174 tells you what the leftmost polynomial's value is at x=0. The second
175 entry discusses the impact of the second polynomial. If you have a
176 vector of p+1 shape function weights, you can compute the scalar
177 product of this vector and this weight vector, and you'll know what
178 the solution at the point x=0 is.
179
180 If the rightmost value, i.e. the value at x=1 is of interest, the order
181 within this vector has to be inverted.
182
183 This is a nodal 1d operation, i.e. there are no integrals or similar
184 involved. We take the real 1d polyomial from the shape function.
185"""
186
187 DerivativeOperator1dDocumentation = """
188 Matrix which is given the p+1 weights of the shape function
189 along a 1d line. It then returns the p+1 gradients at the integration
190 points.
191
192 This is a point-wise operator, i.e. it works directly with the weights
193 of the shape function in the function space. It does not involve any
194 integral et al. Furthermore, it is exact given a certain Lagrange
195 polynomial.
196"""
197
198 InterpolationMatrix1dDocumentation = """
199 This is essentially a @f$ 3 \cdot (order+1)^d-1 @f$ matrix split up
200 into a three-dimensional tensor. The documentation/content of this
201 tensor can be found on the page about DG AMR which is held in exahype
202 within the file exahype2/dg/Riemann.h.
203"""
204
205 RestrictionMatrix1dDocumentation = """
206 This is essentially a @f$ 3 \cdot (order+1)^d-1 @f$ matrix split up
207 into a three-dimensional tensor. The documentation/content of this
208 tensor can be found on the page about DG AMR which is held in exahype
209 within the file exahype2/dg/Riemann.h.
210"""
211
213 """
214
215 See Eq. 3.2 in
216 https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
217
218 """
219 result = [0.0]*self.__num_points
220 for j,xj in enumerate(self.__quadrature_points):
221 divisor = 1.0
222 for k,xk in enumerate(self.__quadrature_points):
223 if j != k:
224 divisor *= (xj-xk)
225 result[j] = 1.0/divisor
226 return result
227
228
230 """
231
232 Reference mass matrix. See class documentation. I moved the docu there,
233 as doxygen does not extract documentation from private routines in Python.
234
235 """
236 mass_matrix = [ [0.0 for _ in range(0,self.dofs_per_axisdofs_per_axis)] for _ in range(0,self.dofs_per_axisdofs_per_axis) ]
237
238 for row in range(0,self.dofs_per_axisdofs_per_axis):
239 for col in range(0,self.dofs_per_axisdofs_per_axis):
240 ans,err = quad(lambda x: self.value1d(x,row) * self.value1d(x,col),0,1)
241 mass_matrix[row][col] = ans
242
243 return [ mass_matrix[i][i] for i in range(0,self.dofs_per_axisdofs_per_axis) ]
244
245
247 """
248
249 Computes the (reference) element stiffness matrix for an approximation of order self.__max_poly_order.
250
251 @see StiffnessOperator1dDocumentation which is also used to dump documentation
252 into generated C++ code.
253
254 """
255 stiffness_matrix = [ [0.0 for _ in range(0,self.dofs_per_axisdofs_per_axis)] for _ in range(0,self.dofs_per_axisdofs_per_axis) ]
256
257 for row in range(0,self.dofs_per_axisdofs_per_axis):
258 for col in range(0,self.dofs_per_axisdofs_per_axis):
259 ans,err = quad(lambda x: self.derivative1d(x,row) * self.value1d(x,col),0,1)
260 stiffness_matrix[row][col] = ans
261
262 return stiffness_matrix
263
264
266 """
267
268 compute the interpolation matrix, which used to prepare halo layer data for fine grid points at the refinement boundary.
269
270 The matrix is only 1d mapping, i.e. used for "surface" of a 2D domain and have a size of 3Dof \times Dof.
271 Its elements are \phi_j(x^f_i), i.e. the value of basis functions for coarse grid on fine grid point.
272 In practive we multiply this matrix to vector of weight in coarse grid (u^c_i):
273
274 \phi_j(x^f_i)*u^c_i=u^c(x^f_i)
275
276 u^c(x) is the original solution on coarse grid. Thus the result is actually the evaluation of coarse solution on fine grid position
277
278 We can also use the same matrix to interpolate along the normal direction.
279
280 """
282 return [[[1.0]],[[1.0]],[[1.0]]]
283 else:
284 interp_matrix = [ [0.0 for _ in range(0,self.dofs_per_axisdofs_per_axis)] for _ in range(0,3*self.dofs_per_axisdofs_per_axis) ]
285 interp_matrix = [ [[0.0 for _ in range(0,self.dofs_per_axisdofs_per_axis)] for _ in range(0,self.dofs_per_axisdofs_per_axis) ] for _ in range(0,3) ]
286 fine_grid_point = [0.0 for _ in range(0,3*self.dofs_per_axisdofs_per_axis) ]
287 for index in range(0,self.dofs_per_axisdofs_per_axis):
288 fine_grid_point[index ]=self.quadrature_points[index]/3.0
289 fine_grid_point[index+ self.dofs_per_axisdofs_per_axis]=self.quadrature_points[index]/3.0+1.0/3.0
290 fine_grid_point[index+2*self.dofs_per_axisdofs_per_axis]=self.quadrature_points[index]/3.0+2.0/3.0
291 for block in range(0, 3):
292 for row in range(0, self.dofs_per_axisdofs_per_axis):
293 for col in range(0, self.dofs_per_axisdofs_per_axis):
294 interp_matrix[block][row][col]=self.value1d(fine_grid_point[row],col)
295
296 return interp_matrix
297
298
300 """
301
302 compute the restriction matrix, which used to prepare halo layer data for coarse cells at the refinement boundary.
303 the matrix formulation depends on the interpolation scheme.
304
305 The matrix is only 1d mapping, i.e. used for "surface" of a 2D domain and have a size of Dof \times 3Dof.
306
307 We can also use the same matrix to interpolate along the normal direction.
308
309
310 """
312 return [[[1.0/3.0]],[[1.0/3.0]],[[1.0/3.0]]]
313 else:
314 restriction_matrix = [ [[0.0 for _ in range(0,self.dofs_per_axisdofs_per_axis)] for _ in range(0,self.dofs_per_axisdofs_per_axis) ] for _ in range(0,3) ]
315 fine_grid_point = [0.0 for _ in range(0,3*self.dofs_per_axisdofs_per_axis) ]
316 for index in range(0,self.dofs_per_axisdofs_per_axis):
317 fine_grid_point[index ]=self.quadrature_points[index]/3.0
318 fine_grid_point[index+ self.dofs_per_axisdofs_per_axis]=self.quadrature_points[index]/3.0+1.0/3.0
319 fine_grid_point[index+2*self.dofs_per_axisdofs_per_axis]=self.quadrature_points[index]/3.0+2.0/3.0
320 for block in range(0, 3):
321 for row in range(0, self.dofs_per_axisdofs_per_axis):
322 for col in range(0, self.dofs_per_axisdofs_per_axis):
323 restriction_matrix[block][row][col]=self.value1d(fine_grid_point[col],row)/3.0
324
325 return restriction_matrix
326
327
328 def __compute_K1(self):
329 """
330
331 @todo Not yet updated atm but seems also not to be used anywhere
332
333 Computes the difference between the reference element mass operator
334 evaluated at point xi=1.0 and the element stiffness matrix.
335
336 :return: delta between the reference element mass operator at point xi=1.0 and the element stiffness matrix
337 """
338 phi1, _ = self.__evaluate(1.0)
339 Kxi = self.__stiffness_operator
340
341 K1 = [[mp.mpf(0) for _ in range(self.__num_points)] for _ in range(self.__num_points)]
342 #FRm = [[mp.mpf(0) for _ in range(self.__num_points)] for _ in range(self.__num_points)]
343
344 for k in range(0, self.__num_points):
345 for l in range(0, self.__num_points):
346 #FRm[k][l] = phi1[k]*phi1[l]
347 K1[k][l] = phi1[k]*phi1[l] - Kxi[k][l]
348
349 #K1_orig = np.subtract(FRm,self.__stiffness_operator)
350 #for k in range(0, self.__num_points):
351 # for l in range(0, self.__num_points):
352 # print(K1[k][l] - K1_orig[k][l])
353
354 return K1
355
356
358 """
359
360 @see DerivativeOperator1dDocumentation which is also used to dump documentation
361 into generated C++ code.
362
363 """
364 dQdx = [ [0.0 for _ in range(0,self.dofs_per_axisdofs_per_axis)] for _ in range(0,self.dofs_per_axisdofs_per_axis) ]
365
366 for row in range(0,self.dofs_per_axisdofs_per_axis):
367 for col in range(0,self.dofs_per_axisdofs_per_axis):
368 dQdx[row][col] = self.derivative1d( self.quadrature_points[row], col)
369
370 return dQdx
371
372
374 """
375
376 @todo Definitely useful, but not used atm
377
378 Transforms the degrees of freedom located on a coarse grid edge
379 quadrature_points to degrees of freedoms located on quadrature_points of a fine grid edge.
380 The difference in levels is 1.
381
382 Let us denote by P the fine grid projector (= equidistantGridProjector). The fine grid DoF
383 are computed according to:
384
385 u^{fine;j}_i = sum_{m} P^{j}_im u^{coarse}_m
386
387 Args:
388 j:
389 subinterval index
390
391 Returns:
392 fineGridProjector:
393 Operator to express polynomial function associated with original interval with basis functions associated with subinterval j
394 """
395 fineGridProjector = [[mp.mpf(0) for _ in range(self.__num_points)] for _ in range(self.__num_points)] # 4 x 4 for self.__max_poly_order=3
396
397 for i in range(0, self.__num_points): # Eq. basis
398 phi_i, _ = self.__evaluate((self.__quadrature_points[i]+j)/mp.mpf(3.0)) # comma after phi_i is important
399 for m in range(0, self.__num_points): # DG basis
400 fineGridProjector[m][i] = phi_i[m]
401 return fineGridProjector
402
403
405 """
406
407 Compute per basis function what the very left value of the polynomial would be
408
409 @see BasisFunctionValuesLeft1dDocumentation which defines what this routine
410 does and is also used to inject documentation into the generated C++ code.
411
412 """
413 return [ self.value1d(0.0, number) for number in range(self.dofs_per_axisdofs_per_axis) ]
414
415
417 """
418
419 @todo Not used atm, but definitely useful
420
421 Transforms the degrees of freedom located at non-equidistant Lagrange support points
422 quadrature_points to degrees of freedoms located at quadrature_points of an equidistant grid over (0,1).
423
424 Let us denote by P the projection operator (= equidistantGridProjector). The equidistant DoF
425 are computed according to:
426
427 u^eq_i = sum_{m} P_im u^DG_m
428
429 Returns:
430 equidistantGridProjector:
431 The corresponding degrees of freedom located at quadrature_points of an equidistant grid over (0,1).
432 """
433 equidistantGridProjector = [[mp.mpf(0) for _ in range(self.__num_points)] for _ in range(self.__num_points)] # 4 x 4 for self.__max_poly_order=3
434 subxi = mp.linspace(mp.mpf(0.0), mp.mpf(1.0), self.__num_points)
435
436 for i in range(0, self.__num_points): # Eq. basis
437 phi_i, _ = self.__evaluate(subxi[i]) # comma after phi_i is important
438 for m in range(0, self.__num_points): # DG basis
439 equidistantGridProjector[m][i] = phi_i[m]
440 return equidistantGridProjector
441
442
443 def init_dictionary_with_default_parameters(self,dictionary,use_multidimensional_arrays):
444 def snake_to_camel(word):
445 return ''.join(x.capitalize() or '_' for x in word.lower().split('_'))
446
447 basisDeclarations = ""
448 basisInitializers = ""
449
450 dictionary["ORDER"] = self.order
451 dictionary["BASIS_DECLARATIONS"] = ""
452 dictionary["BASIS_INITIALIZERS"] = ""
453
454 dictionary["BASIS_DECLARATIONS"] += "/**\n"
455 dictionary["BASIS_DECLARATIONS"] += self.QuadraturePoints1dDocumentation
456 dictionary["BASIS_DECLARATIONS"] += "*/\n"
457 #dictionary["BASIS_DECLARATIONS"] += "static const double QuadraturePoints1d[DGOrder+1];\n\n"
458 #dictionary["BASIS_INITIALIZERS"] += " QuadraturePoints1d{},\n".format( render_tensor_1(self.quadrature_points) )
459 dictionary["BASIS_DECLARATIONS"] += "static constexpr double QuadraturePoints1d[] = {};\n\n".format( render_tensor_1(self.quadrature_points) )
460
461 dictionary["BASIS_DECLARATIONS"] += "/**\n"
462 dictionary["BASIS_DECLARATIONS"] += self.QuadratureWeights1dDocumentation
463 dictionary["BASIS_DECLARATIONS"] += "*/\n"
464 #dictionary["BASIS_DECLARATIONS"] += "const double QuadratureWeights1d[DGOrder+1];\n\n"
465 #dictionary["BASIS_INITIALIZERS"] += " QuadratureWeights1d{},\n".format( render_tensor_1(self.quadrature_weights) )
466 dictionary["BASIS_DECLARATIONS"] += "static constexpr double QuadratureWeights1d[] = {};\n\n".format( render_tensor_1(self.quadrature_weights) )
467
468 dictionary["BASIS_DECLARATIONS"] += "/**\n"
469 dictionary["BASIS_DECLARATIONS"] += self.MassMatrix1dDocumentation
470 dictionary["BASIS_DECLARATIONS"] += "*/\n"
471 #dictionary["BASIS_DECLARATIONS"] += "const double MassMatrixDiagonal1d[DGOrder+1];\n\n"
472 #dictionary["BASIS_INITIALIZERS"] += " MassMatrixDiagonal1d{},\n".format( render_tensor_1(self.__compute_mass_matrix_diagonal()) )
473 dictionary["BASIS_DECLARATIONS"] += "static constexpr double MassMatrixDiagonal1d[] = {};\n\n".format( render_tensor_1(self.__compute_mass_matrix_diagonal()) )
474
475 dictionary["BASIS_DECLARATIONS"] += "/**\n"
476 dictionary["BASIS_DECLARATIONS"] += self.StiffnessOperator1dDocumentation
477 dictionary["BASIS_DECLARATIONS"] += "*/\n"
478 if use_multidimensional_arrays:
479 #dictionary["BASIS_DECLARATIONS"] += "const double StiffnessOperator1d[DGOrder+1][DGOrder+1];\n\n"
480 #dictionary["BASIS_INITIALIZERS"] += " StiffnessOperator1d{},\n".format( render_tensor_2(self.__compute_stiffness_operator(),use_multidimensional_arrays) )
481 dictionary["BASIS_DECLARATIONS"] += "static constexpr double StiffnessOperator1d[][] = {};\n\n".format( render_tensor_2(self.__compute_stiffness_operator(),use_multidimensional_arrays) )
482 else:
483 #dictionary["BASIS_DECLARATIONS"] += "const double StiffnessOperator1d[(DGOrder+1)*(DGOrder+1)];\n\n"
484 #dictionary["BASIS_INITIALIZERS"] += " StiffnessOperator1d{},\n".format( render_tensor_2(self.__compute_stiffness_operator(),use_multidimensional_arrays) )
485 dictionary["BASIS_DECLARATIONS"] += "static constexpr double StiffnessOperator1d[] = {};\n\n".format( render_tensor_2(self.__compute_stiffness_operator(),use_multidimensional_arrays) )
486
487 dictionary["BASIS_DECLARATIONS"] += "/**\n"
488 dictionary["BASIS_DECLARATIONS"] += self.BasisFunctionValuesLeft1dDocumentation
489 dictionary["BASIS_DECLARATIONS"] += "*/\n"
490 dictionary["BASIS_DECLARATIONS"] += "static constexpr double BasisFunctionValuesLeft1d[] = {};\n\n".format( render_tensor_1(self.__compute_basis_function_values_left()) )
491
492 dictionary["BASIS_DECLARATIONS"] += "/**\n"
493 dictionary["BASIS_DECLARATIONS"] += self.DerivativeOperator1dDocumentation
494 dictionary["BASIS_DECLARATIONS"] += "*/\n"
495 if use_multidimensional_arrays:
496 dictionary["BASIS_DECLARATIONS"] += "static constexpr double DerivativeOperator1d[][] = {};\n\n".format( render_tensor_2(self.__compute_derivative_operator(),use_multidimensional_arrays) )
497 else:
498 dictionary["BASIS_DECLARATIONS"] += "static constexpr double DerivativeOperator1d[] = {};\n\n".format( render_tensor_2(self.__compute_derivative_operator(),use_multidimensional_arrays) )
499
500 dictionary["BASIS_DECLARATIONS"] += "/**\n"
501 dictionary["BASIS_DECLARATIONS"] += self.RestrictionMatrix1dDocumentation
502 dictionary["BASIS_DECLARATIONS"] += "*/\n"
503 if use_multidimensional_arrays:
504 dictionary["BASIS_DECLARATIONS"] += "static constexpr double RestrictionMatrix1d[][][] = {};\n\n".format( render_tensor_3(self.__compute_restriction_matrix(),use_multidimensional_arrays) )
505 else:
506 dictionary["BASIS_DECLARATIONS"] += "static constexpr double RestrictionMatrix1d[] = {};\n\n".format( render_tensor_3(self.__compute_restriction_matrix(),use_multidimensional_arrays) )
507
508 dictionary["BASIS_DECLARATIONS"] += "/**\n"
509 dictionary["BASIS_DECLARATIONS"] += self.InterpolationMatrix1dDocumentation
510 dictionary["BASIS_DECLARATIONS"] += "*/\n"
511 if use_multidimensional_arrays:
512 dictionary["BASIS_DECLARATIONS"] += "static constexpr double InterpolationMatrix1d[][][] = {};\n\n".format( render_tensor_3(self.__compute_interpolation_matrix(),use_multidimensional_arrays) )
513 else:
514 dictionary["BASIS_DECLARATIONS"] += "static constexpr double InterpolationMatrix1d[] = {};\n\n".format( render_tensor_3(self.__compute_interpolation_matrix(),use_multidimensional_arrays) )
515
516
517 @property
518 @abstractmethod
520 assert False, "to be implemented by subclass"
521
522
523 @property
524 @abstractmethod
526 assert False, "to be implemented by subclass"
527
528
529 def value1d(self,x,number):
530 """
531
532 Return the numberth shape function's value at x
533
534 This equals Dominic's old evaluate function, but that one evaluated all the
535 polynomials in one rush, whereas this one only looks at shape function number
536 number. Also, Dominic's version did both the value and the derivative. I
537 split this one up into two version.
538
539 """
540 enumerator = 1.0
541 for i in range(0,self.dofs_per_axisdofs_per_axis):
542 if i!=number:
543 enumerator *= (x-self.quadrature_points[i])
544
545 denominator = 1.0
546 for i in range(0,self.dofs_per_axisdofs_per_axis):
547 if i!=number:
548 denominator *= (self.quadrature_points[number]-self.quadrature_points[i])
549
550 return enumerator/denominator
551
552
553 def derivative1d(self,x,number):
554 """
555
556 Return the numberth shape function's derivative at x
557
558 u = sum_i u_i phi_i
559
560 => d/dx u = sum_j (d/dx phi_j) u_j
561
562 To construct the gradient field, we make the ansatz:
563
564 grad u = sum_i g_i phi_i
565
566 =>
567
568 g_i = ( d/dx u, phi_i )_[0,1] / (phi_i,phi_i)_[0,1] = sum_j (d/dx phi_j, phi_i)_[0,1] / w_i
569
570 = w_i (d/dx phi_j) (x_i) / w_i = (d/dx phi_j) (x_i)
571
572 = DUDX^T
573
574 where DUDX is the operator computed by this function:
575
576 DUDX_ij = (d/dx phi_i) (x_j)
577
578 It can be further written as
579
580 DUDX_ij = 1/w_i * K^T_ij
581
582 where the stiffness matrix K is defined as
583
584 K_ij = <d/dx phi_i, phi_j>_[0,1] = w_j (d/dx phi_j) (x_i)
585
586 :return: transposed derivative operator
587
588 :note: If you want to use this operator to compute the gradient of the solution,
589 you need to use the transpose.
590
591
592 """
593
594 result = 0.0
595 for l in range(0,self.dofs_per_axisdofs_per_axis):
596 if l!=number:
597 enumerator = 1.0
598 denominator = 1.0/(self.quadrature_points[number]-self.quadrature_points[l])
599 for m in range(0,self.dofs_per_axisdofs_per_axis):
600 if m!=number and m!=l:
601 enumerator *= x-self.quadrature_points[m]
602 denominator *= 1.0/(self.quadrature_points[number]-self.quadrature_points[m])
603 result += enumerator*denominator
604 return result
605
606from scipy.special import legendre
607from scipy.special import roots_sh_legendre
608
610 """
611
612 The Gauss-Legendre Basis is by construction the only basis which
613 yields diagonal mass matrices. I rely solely on SciPy here to
614 construct these functions, as I'm interested in numerical values
615 only, and as the normal SciPy accuracy is sufficient.
616
617 There might be more elegant implementations using SymPy which
618 exploit the special structure behind Gauss-Legendre. But I prefer
619 here to stick to one Python package only.
620
621 """
622
623 def __init__(self,polynomial_order):
624 super(GaussLegendreBasis,self).__init__(polynomial_order)
625
626
627 @property
629 """
630
631 It seems to be a little bit weird that I have to use the polynomial
632 order plus one here, but you have to read the docu to understand
633 that we are evaluating the polynomial of a certain order to get the
634 roots. So the argument has to be equal to the number of quadrature
635 points that I wannt have.
636
637 I originally used the roots_legendre and then did recalibrate the
638 result manually. Today, SciPy directly offers a shifted Legendre
639 function which I can use directly. So I don't need an expression like
640
641
642 return [ x/2.0+0.5 for x in roots_on_minus_one_one_interval ]
643
644
645 anymore. The nicest definition of the shifted Legendre shape functions
646 can, as so often, be found on Wikipedia:
647
648 https://en.wikipedia.org/wiki/Legendre_polynomials#Shifted_Legendre_polynomials
649
650 """
651 roots, weights = scipy.special.roots_sh_legendre(self._polynomial_order_polynomial_order+1)
652 return roots
653
654
655 @property
657 """
658
659 It seems to be a little bit weird that I have to use the polynomial
660 order plus one here, but you have to read the docu to understand
661 that we are evaluating the polynomial of a certain order to get the
662 roots. So the argument has to be equal to the number of quadrature
663 points that I wannt have.
664
665 I originally used the roots_legendre and then did recalibrate the
666 result manually. Today, SciPy directly offers a shifted Legendre
667 function which I can use directly. So I don't need an expression like
668
669
670 return [ x/2.0+0.5 for x in roots_on_minus_one_one_interval ]
671
672
673 anymore. The nicest definition of the shifted Legendre shape functions
674 can, as so often, be found on Wikipedia:
675
676 https://en.wikipedia.org/wiki/Legendre_polynomials#Shifted_Legendre_polynomials
677
678 """
679 roots, weights = scipy.special.roots_sh_legendre(self._polynomial_order_polynomial_order+1)
680 return weights
681
682
683
685 """
686
687 The Gauss-Lobatto quadrature points do not yield a diagonal mass
688 matrix. They yield a full matrix. However, there are codes which
689 ignore this fact and lump the basis. Mathematically, this is
690 equivalent to the Gauss-Lobatto quadrature because this leads to a diagonal matrix that is easy to invert
691
692 """
693
694 def __init__(self,polynomial_order):
695 super(GaussLobattoBasisWithLumpedDiagonalBasis,self).__init__(polynomial_order)
696
697
698 @property
700 assert False, "to be implemented by subclass"
701 raise Exception( "check implementation. See class comment")
702 x,w = quadrature.gauss_lobatto(self.dofs_per_axisdofs_per_axis,PREC_DIGITS)
703 return self._transform_and_sort_points_and_weights(x,w)
704
The Gauss-Legendre Basis is by construction the only basis which yields diagonal mass matrices.
quadrature_weights(self)
It seems to be a little bit weird that I have to use the polynomial order plus one here,...
quadrature_points(self)
It seems to be a little bit weird that I have to use the polynomial order plus one here,...
__compute_stiffness_operator(self)
Computes the (reference) element stiffness matrix for an approximation of order self....
__compute_restriction_matrix(self)
compute the restriction matrix, which used to prepare halo layer data for coarse cells at the refinem...
__compute_interpolation_matrix(self)
compute the interpolation matrix, which used to prepare halo layer data for fine grid points at the r...
__compute_basis_function_values_left(self)
Compute per basis function what the very left value of the polynomial would be.
init_dictionary_with_default_parameters(self, dictionary, use_multidimensional_arrays)
To be implemented by subclass.
Abstract base class for all Lagrange bases available in ExaHyPE 2.
init_dictionary_with_default_parameters(self, dictionary, use_multidimensional_arrays)
To be implemented by subclass.