diff --git a/.gitignore b/.gitignore index e00918d..d8beb52 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ __pycache__ -*.lock +# don't commit the lock; it may get outdated to quickly +poetry.lock + diff --git a/README.md b/README.md index dbb5f5c..4240815 100644 --- a/README.md +++ b/README.md @@ -23,3 +23,16 @@ If we use basis functions for order $\mathcal{O}$, we have $\frac{1}{6} \times \ | $M_{kl}$ | $\int_T \Psi_k \Psi_l dx$ | `M3` | | | $\int_T \Phi_k \Phi_l dx$ | `M2` | | $F_{kl}^{-,j}$ | $\int_T \Psi_k \Psi_l dx$ | `rT` | + +## Usage + +Use poetry to install or use the package. In particular, consider the commands + +``` +# install +poetry install +# format files +poetry run black src +# run unit tests +poetry run coverage run -m nose2 +``` diff --git a/src/seissol_matrices/base.py b/src/seissol_matrices/base.py new file mode 100644 index 0000000..2137683 --- /dev/null +++ b/src/seissol_matrices/base.py @@ -0,0 +1,22 @@ +import numpy as np + + +def collocate(basis, points): + # takes a basis object, and a points array + # of the form: npoints × dim + # outputs a collocation matrix of the form + # npoints × nbasis + + # as far as the author of these lines is aware, + # this module has no broadcasing functionality yet + # for basis functions. Which is sad. + + assert basis.dim() == points.shape[1] + + nbasis = basis.number_of_basis_functions() + + coll = np.empty((nbasis, points.shape[0])) + for i in range(nbasis): + for j in range(points.shape[0]): + coll[i, j] = basis.eval_basis(points[j, :], i) + return coll diff --git a/src/seissol_matrices/basis_functions.py b/src/seissol_matrices/basis_functions.py index 9c6a107..402d3c8 100644 --- a/src/seissol_matrices/basis_functions.py +++ b/src/seissol_matrices/basis_functions.py @@ -40,7 +40,7 @@ def singularity_free_jacobi_polynomial(n, a, b, x, y): def singularity_free_jacobi_polynomial_and_derivatives(n, a, b, x, y): if n == 0: - return 1.0, 0.0, 0.0 + return np.ones(np.shape(x)), np.zeros(np.shape(x)), np.zeros(np.shape(x)) pm_1 = 1.0 ddx_pm_1 = 0.0 ddy_pm_1 = 0.0 @@ -85,11 +85,17 @@ def eval_diff_basis(self, x, i, k): def number_of_basis_functions(self, o): pass + def dim(self): + pass + ################################################################################ class BasisFunctionGenerator1D(BasisFunctionGenerator): + def dim(self): + return 1 + def eval_basis(self, x, i): r_num = 2 * x - 1.0 r_den = 1.0 @@ -97,7 +103,7 @@ def eval_basis(self, x, i): def eval_diff_basis(self, x, i, k): if i == 0: - return 0 + return np.zeros(x.shape) else: r_num = 2 * x - 1.0 r_den = 1.0 @@ -113,6 +119,9 @@ def number_of_basis_functions(self): class BasisFunctionGenerator2D(BasisFunctionGenerator): + def dim(self): + return 2 + def unroll_index(self, i): n = i[0] + i[1] tri = 0.5 * n * (n + 1) @@ -171,6 +180,9 @@ def number_of_basis_functions(self): class BasisFunctionGenerator3D(BasisFunctionGenerator): + def dim(self): + return 3 + def unroll_index(self, i): n = i[0] + i[1] + i[2] tet = (n * (n + 1) * (n + 2)) / 6.0 diff --git a/src/seissol_matrices/dg_matrices.py b/src/seissol_matrices/dg_matrices.py index 05c71f1..114d6a2 100644 --- a/src/seissol_matrices/dg_matrices.py +++ b/src/seissol_matrices/dg_matrices.py @@ -5,6 +5,10 @@ import quad_rules.GaussJacobi import quad_rules.Quadrature +import itertools + +from seissol_matrices import base + from seissol_matrices import basis_functions @@ -16,6 +20,9 @@ def __init__(self, o, d): self.order = o self.dim = d + self.M = None + self.K = self.dim * [None] + if self.dim == 3: self.generator = basis_functions.BasisFunctionGenerator3D(self.order) self.geometry = np.array( @@ -27,71 +34,158 @@ def __init__(self, o, d): ] ) self.face_generator = dg_generator(o, 2) - n, w = quad_rules.JaskowiecSukumar.JaskowiecSukumar().find_best_rule( - 2 * self.order + self.quadrule_finder = ( + quad_rules.JaskowiecSukumar.JaskowiecSukumar().find_best_rule ) + self.generator_finder = basis_functions.BasisFunctionGenerator3D elif self.dim == 2: self.generator = basis_functions.BasisFunctionGenerator2D(self.order) self.geometry = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]) - n, w = quad_rules.WitherdenVincentTri.WitherdenVincentTri().find_best_rule( - 2 * self.order + self.quadrule_finder = ( + quad_rules.WitherdenVincentTri.WitherdenVincentTri().find_best_rule ) + self.generator_finder = basis_functions.BasisFunctionGenerator2D elif self.dim == 1: self.generator = basis_functions.BasisFunctionGenerator1D(self.order) self.geometry = np.array([[0.0], [1.0]]) - n, w = quad_rules.GaussJacobi.GaussJacobi(0, 0).find_best_rule( - 2 * self.order - ) + self.quadrule_finder = quad_rules.GaussJacobi.GaussJacobi( + 0, 0 + ).find_best_rule + self.generator_finder = basis_functions.BasisFunctionGenerator1D else: raise Execption("Can only generate 1D, 2D or 2D basis functions") - self.nodes, self.weights = quad_rules.Quadrature.transform(n, w, self.geometry) - self.M = None - self.K = self.dim * [None] + self.nodes, self.weights = self.get_quadrature_rule(2 * self.order) + + def get_quadrature_rule(self, order): + n, w = self.quadrule_finder(order) + + return quad_rules.Quadrature.transform(n, w, self.geometry) + + def multilinear_form(self, face, der, order=None, dim=None, side=None): + """ + Given an array der, we compute a tensor V[i...] of dimension + len(der) defined by (given in C++ parameter pack notation) + + V[i...] = integral(1 * ... * f_i^(der)) + + For example, to obtain a mass matrix, choose der = [None,None]. + For a stiffness matrix in x-direction, choose der = [0,None]. + For a trilinear form, choose der = [None,None,None]. + + Optionally can supply different orders for the basis functions. (using the order parameter) + If not given, the order will assumed to be self.order for all constituents. + + Moreover, we can opt to take a projected face basis instead of a volume basis via the + side argument. If an entry is a negative number there, we assume a non-rotated face, + otherwise the respective face side. If an entry is None, then we assume the respective volume term. + """ + + # NOTE: we could extract the multiplication/quadrature part from the function selection part. + # I.e., eval_functions (or a list of basis functions here) could be given as an input parameter, + # instead of the two current ones. + + # we only support a first derivative at the moment + assert np.all([d in (None, *(i for i in range(self.dim))) for d in der]) + + if order is None: + order = [self.order] * len(der) + + if side is None: + side = [None] * len(der) + + if dim is None: + dim = [2] * len(der) if face else [3] * len(der) + + assert len(order) == len(der) + assert len(side) == len(der) + + def get_basis_function(o, d, s, r): + if r == 3: + basis = self.generator_finder(o) + else: + basis = self.face_generator.generator_finder(o) + basiseval_pre = ( + basis.eval_basis + if d is None + else lambda x, i: basis.eval_diff_basis(x, i, d) + ) + if face: + if r == 3: + basiseval = lambda x, i: basiseval_pre( + self.volume_to_face_parametrisation(x, s), i + ) + else: + if s is None: + basiseval = basiseval_pre + else: + basiseval = lambda x, i: basiseval_pre( + self.face_to_face_parametrisation(x, s), i + ) + else: + # volume * face does not make sense when evaluating on the whole volume + assert r == 3 + basiseval = basiseval_pre + return basiseval, basis.number_of_basis_functions() + + basis_functions = [ + get_basis_function(o, d, s, r) for o, d, s, r in zip(order, der, side, dim) + ] + + if face: + nodes, weights = self.face_generator.get_quadrature_rule(np.sum(order)) + else: + nodes, weights = self.get_quadrature_rule(np.sum(order)) + + sizes = [bn for _, bn in basis_functions] + tensor = np.empty(sizes) + generic_eval_basis = lambda x, index: np.prod( + [basis_functions[k][0](x, i) for k, i in enumerate(index)], axis=0 + ) + for index in itertools.product(*[range(size) for size in sizes]): + eval_basis = lambda x: generic_eval_basis(x, index) + tensor[index] = quad_rules.Quadrature.quad(nodes, weights, eval_basis) + return tensor def mass_matrix(self): if not np.any(self.M == None): return self.M - number_of_basis_functions = self.generator.number_of_basis_functions() - self.M = np.zeros((number_of_basis_functions, number_of_basis_functions)) - - for i in range(number_of_basis_functions): - for j in range(number_of_basis_functions): - prod = lambda x: self.generator.eval_basis( - x, i - ) * self.generator.eval_basis(x, j) - self.M[i, j] = quad_rules.Quadrature.quad( - self.nodes, self.weights, prod - ) + + self.M = self.multilinear_form(False, [None, None]) return self.M def stiffness_matrix(self, dim): if not np.any(self.K[dim] == None): return self.K[dim] - number_of_basis_functions = self.generator.number_of_basis_functions() - self.K[dim] = np.zeros((number_of_basis_functions, number_of_basis_functions)) - - for i in range(number_of_basis_functions): - for j in range(number_of_basis_functions): - prod = lambda x: self.generator.eval_diff_basis( - x, i, dim - ) * self.generator.eval_basis(x, j) - self.K[dim][i, j] = quad_rules.Quadrature.quad( - self.nodes, self.weights, prod - ) + + self.K[dim] = self.multilinear_form(False, [dim, None]) return self.K[dim] - def kDivM(self, dim): - stiffness = self.stiffness_matrix(dim) + def kDivM(self, dim, matorder=None): mass = self.mass_matrix() - + if matorder is None: + stiffness = self.stiffness_matrix(dim) + else: + stiffness = self.multilinear_form( + False, [None, dim, None], order=[matorder, self.order, self.order] + ) return np.linalg.solve(mass, stiffness) - def kDivMT(self, dim): - stiffness = self.stiffness_matrix(dim) + def kDivMT(self, dim, matorder=None): mass = self.mass_matrix() + if matorder is None: + stiffness = self.stiffness_matrix(dim) + sT = stiffness.T + else: + correct1 = self.multilinear_form( + False, [dim, None, None], order=[matorder, self.order, self.order] + ) + correct2 = self.multilinear_form( + False, [None, None, dim], order=[matorder, self.order, self.order] + ) + sT = correct1 + correct2 - return np.linalg.solve(mass, stiffness.T) + return np.linalg.solve(mass, sT) def face_to_face_parametrisation(self, x, side): # implement Dumbser & Käser, 2006 Table 2 b) @@ -113,28 +207,7 @@ def face_to_face_parametrisation(self, x, side): def face_times_face_mass_matrix(self, side): assert self.dim == 3 - number_of_face_basis_functions = ( - self.face_generator.generator.number_of_basis_functions() - ) - - projected_basis_function = ( - lambda x, i: self.face_generator.generator.eval_basis( - self.face_to_face_parametrisation(x, side), i - ) - ) - - matrix = np.zeros( - (number_of_face_basis_functions, number_of_face_basis_functions) - ) - for i in range(number_of_face_basis_functions): - for j in range(number_of_face_basis_functions): - prod = lambda x: self.face_generator.generator.eval_basis( - x, i - ) * projected_basis_function(x, j) - matrix[i, j] = quad_rules.Quadrature.quad( - self.face_generator.nodes, self.face_generator.weights, prod - ) - return matrix + return self.multilinear_form(True, [None, None], side=[None, side], dim=[2, 2]) def fP(self, side): return self.face_times_face_mass_matrix(side) @@ -164,25 +237,7 @@ def volume_to_face_parametrisation(self, x, side): def volume_times_face_mass_matrix(self, side): assert self.dim == 3 - number_of_face_basis_functions = ( - self.face_generator.generator.number_of_basis_functions() - ) - - projected_basis_function = lambda x, i: self.generator.eval_basis( - self.volume_to_face_parametrisation(x, side), i - ) - number_of_basis_functions = self.generator.number_of_basis_functions() - - matrix = np.zeros((number_of_face_basis_functions, number_of_basis_functions)) - for i in range(number_of_face_basis_functions): - for j in range(number_of_basis_functions): - prod = lambda x: self.face_generator.generator.eval_basis( - x, i - ) * projected_basis_function(x, j) - matrix[i, j] = quad_rules.Quadrature.quad( - self.face_generator.nodes, self.face_generator.weights, prod - ) - return matrix + return self.multilinear_form(True, [None, None], side=[None, side], dim=[2, 3]) def fMrT(self, side): return self.volume_times_face_mass_matrix(side) @@ -192,11 +247,32 @@ def rT(self, side): mass = self.face_generator.mass_matrix() return np.linalg.solve(mass, matrix) - def rDivM(self, side): + def rDivM(self, side, matorder=None): assert self.dim == 3 - matrix = self.rT(side) mass = self.mass_matrix() - return np.linalg.solve(mass, matrix.T) + if matorder is None: + matrix = self.rT(side).T + else: + prematrix = self.multilinear_form( + True, + [None, None, None], + order=[matorder, self.order, self.order], + side=[None, None, side], + dim=[3, 2, 3], + ) + facemass = self.face_generator.mass_matrix() + matrix = np.linalg.solve(facemass, prematrix).transpose((0, 2, 1)) + return np.linalg.solve(mass, matrix) + + def collocate_volume(self, points): + return base.collocate(self.generator, points) + + def collocate_face(self, points, side): + # points are meant to be 2D here + + # this method wants dim × npoints; but points is given the other way. So, we transpose it twice + projected = self.volume_to_face_parametrisation(points.T, side).T + return base.collocate(self.generator, projected) if __name__ == "__main__": diff --git a/src/seissol_matrices/dr_matrices.py b/src/seissol_matrices/dr_matrices.py index be1b7c8..a22ed1b 100755 --- a/src/seissol_matrices/dr_matrices.py +++ b/src/seissol_matrices/dr_matrices.py @@ -15,70 +15,63 @@ def __init__(self, o, q): self.dg3_generator = dg_matrices.dg_generator(self.order, 3) self.dg2_generator = dg_matrices.dg_generator(self.order, 2) - def V3mTo2n(self, a, b): + def V3mTo2n(self, a, b, matorder=None): m = self.bf3_generator.number_of_basis_functions() points = self.quadrule.points() n = points.shape[0] - matrix = np.zeros((n, m)) + if matorder is None: + evalfun = lambda q, j, k: self.bf3_generator.eval_basis(q, j) + materialdim = 1 + else: + materialbasis = basis_functions.BasisFunctionGenerator3D(matorder) + evalfun = lambda q, j, k: self.bf3_generator.eval_basis( + q, j + ) * materialbasis.eval_basis(q, k) + materialdim = materialbasis.number_of_basis_functions() + + matrix = np.zeros((materialdim, n, m)) for i, p in enumerate(points): for j in range(m): - if a == 0 and b == 0: - matrix[i, j] = self.bf3_generator.eval_basis([p[1], p[0], 0], j) - if a == 0 and b == 1: - matrix[i, j] = self.bf3_generator.eval_basis([p[0], p[1], 0], j) - if a == 0 and b == 2: - matrix[i, j] = self.bf3_generator.eval_basis( - [p[1], 1 - p[0] - p[1], 0], j - ) - if a == 0 and b == 3: - matrix[i, j] = self.bf3_generator.eval_basis( - [1 - p[0] - p[1], p[0], 0], j - ) - if a == 1 and b == 0: - matrix[i, j] = self.bf3_generator.eval_basis([p[0], 0, p[1]], j) - if a == 1 and b == 1: - matrix[i, j] = self.bf3_generator.eval_basis([p[1], 0, p[0]], j) - if a == 1 and b == 2: - matrix[i, j] = self.bf3_generator.eval_basis( - [1 - p[0] - p[1], 0, p[1]], j - ) - if a == 1 and b == 3: - matrix[i, j] = self.bf3_generator.eval_basis( - [p[0], 0, 1 - p[0] - p[1]], j - ) - if a == 2 and b == 0: - matrix[i, j] = self.bf3_generator.eval_basis([0, p[1], p[0]], j) - if a == 2 and b == 1: - matrix[i, j] = self.bf3_generator.eval_basis([0, p[0], p[1]], j) - if a == 2 and b == 2: - matrix[i, j] = self.bf3_generator.eval_basis( - [0, p[1], 1 - p[0] - p[1]], j - ) - if a == 2 and b == 3: - matrix[i, j] = self.bf3_generator.eval_basis( - [0, 1 - p[0] - p[1], p[0]], j - ) - if a == 3 and b == 0: - matrix[i, j] = self.bf3_generator.eval_basis( - [1 - p[0] - p[1], p[0], p[1]], j - ) - if a == 3 and b == 1: - matrix[i, j] = self.bf3_generator.eval_basis( - [1 - p[0] - p[1], p[1], p[0]], j - ) - if a == 3 and b == 2: - matrix[i, j] = self.bf3_generator.eval_basis( - [p[0], 1 - p[0] - p[1], p[1]], j - ) - if a == 3 and b == 3: - matrix[i, j] = self.bf3_generator.eval_basis( - [p[1], p[0], 1 - p[0] - p[1]], j - ) + for k in range(materialdim): + if a == 0 and b == 0: + matrix[k, i, j] = evalfun([p[1], p[0], 0], j, k) + if a == 0 and b == 1: + matrix[k, i, j] = evalfun([p[0], p[1], 0], j, k) + if a == 0 and b == 2: + matrix[k, i, j] = evalfun([p[1], 1 - p[0] - p[1], 0], j, k) + if a == 0 and b == 3: + matrix[k, i, j] = evalfun([1 - p[0] - p[1], p[0], 0], j, k) + if a == 1 and b == 0: + matrix[k, i, j] = evalfun([p[0], 0, p[1]], j, k) + if a == 1 and b == 1: + matrix[k, i, j] = evalfun([p[1], 0, p[0]], j, k) + if a == 1 and b == 2: + matrix[k, i, j] = evalfun([1 - p[0] - p[1], 0, p[1]], j, k) + if a == 1 and b == 3: + matrix[k, i, j] = evalfun([p[0], 0, 1 - p[0] - p[1]], j, k) + if a == 2 and b == 0: + matrix[k, i, j] = evalfun([0, p[1], p[0]], j, k) + if a == 2 and b == 1: + matrix[k, i, j] = evalfun([0, p[0], p[1]], j, k) + if a == 2 and b == 2: + matrix[k, i, j] = evalfun([0, p[1], 1 - p[0] - p[1]], j, k) + if a == 2 and b == 3: + matrix[k, i, j] = evalfun([0, 1 - p[0] - p[1], p[0]], j, k) + if a == 3 and b == 0: + matrix[k, i, j] = evalfun([1 - p[0] - p[1], p[0], p[1]], j, k) + if a == 3 and b == 1: + matrix[k, i, j] = evalfun([1 - p[0] - p[1], p[1], p[0]], j, k) + if a == 3 and b == 2: + matrix[k, i, j] = evalfun([p[0], 1 - p[0] - p[1], p[1]], j, k) + if a == 3 and b == 3: + matrix[k, i, j] = evalfun([p[1], p[0], 1 - p[0] - p[1]], j, k) + if matorder is None: + matrix = np.squeeze(matrix) return matrix - def V3mTo2nTWDivM(self, a, b): - matrix = self.V3mTo2n(a, b) + def V3mTo2nTWDivM(self, a, b, matorder=None): + matrix = self.V3mTo2n(a, b, matorder) mass = self.dg3_generator.mass_matrix() weights = self.quadrule.weights() n = weights.shape[0] @@ -86,7 +79,9 @@ def V3mTo2nTWDivM(self, a, b): for i in range(n): W[i, i] = weights[i] - return np.linalg.solve(mass, np.dot(matrix.T, W)) + matrixT = matrix.T if matorder is None else matrix.transpose((0, 2, 1)) + + return np.linalg.solve(mass, np.dot(matrixT, W)) def quadpoints(self): points = self.quadrule.points() diff --git a/src/seissol_matrices/generate.py b/src/seissol_matrices/generate.py new file mode 100644 index 0000000..5945204 --- /dev/null +++ b/src/seissol_matrices/generate.py @@ -0,0 +1,142 @@ +import sys +from seissol_matrices import json_io +from dg_matrices import dg_generator +from dr_matrices import dr_generator +from vtk_points import vtk_lagrange_2d, vtk_lagrange_3d +from seissol_matrices import quad_points +from plasticity_matrices import PlasticityGenerator + + +def main(): + assert len(sys.argv) > 1 + if sys.argv[1] == "vtk": + import dg_matrices + + vtkpoints2 = {} + vtkpoints3 = {} + for deg in range(0, 9): + vtkpoints2[deg] = vtk_lagrange_2d(deg) + vtkpoints3[deg] = vtk_lagrange_3d(deg) + json_io.write_tensor(vtkpoints2[deg].T, f"vtk2d({deg})", f"vtkbase.json") + json_io.write_tensor(vtkpoints3[deg].T, f"vtk3d({deg})", f"vtkbase.json") + for basisorder in range(2, 9): + dggen = dg_generator(basisorder, 3) + for deg in range(0, 9): + json_io.write_tensor( + dggen.collocate_volume(vtkpoints3[deg]).T, + f"collvv({basisorder},{deg})", + f"vtko{basisorder}.json", + ) + json_io.write_tensor( + dggen.face_generator.collocate_volume(vtkpoints2[deg]).T, + f"collff({basisorder},{deg})", + f"vtko{basisorder}.json", + ) + for f in range(4): + json_io.write_tensor( + dggen.collocate_face(vtkpoints2[deg], f).T, + f"collvf({basisorder},{deg},{f})", + f"vtko{basisorder}.json", + ) + elif sys.argv[1] == "ader": + materialorder = None + if len(sys.argv) > 2: + materialorder = int(sys.argv[2]) + for order in range(2, 9): + filename = f"linear-ader-{order}-h{materialorder}.json" + if materialorder is None: + filename = f"linear-ader-{order}.json" + dggen = dg_generator(order, 3) + for dim in range(3): + json_io.write_tensor( + dggen.kDivMT(dim, materialorder), + f"kDivMT({dim})", + filename, + ) + json_io.write_tensor( + dggen.kDivM(dim, materialorder), f"kDivM({dim})", filename + ) + for line in range(3): + json_io.write_tensor(dggen.fP(line), f"fP({line})", filename) + for face in range(4): + json_io.write_tensor( + dggen.rDivM(face, materialorder), f"rDivM({face})", filename + ) + json_io.write_tensor(dggen.rT(face), f"rT({face})", filename) + json_io.write_tensor(dggen.fMrT(face), f"fMrT({face})", filename) + elif sys.argv[1] == "dr": + materialorder = None + if len(sys.argv) > 2: + materialorder = int(sys.argv[2]) + for order in range(2, 9): + for quadname, quadrule in zip( + ["stroud", "dunavant", "witherden_vincent"], + [ + quad_points.stroud(order + 1), + quad_points.dunavant(order + 1), + quad_points.witherden_vincent(order + 1), + ], + ): + filename = f"dr-{quadname}-{order}-h{materialorder}.json" + if materialorder is None: + filename = f"dr-{quadname}-{order}.json" + generator = dr_generator(order, quadrule) + + quadpoints = generator.quadpoints() + quadweights = generator.quadweights() + resample = generator.resample() + V2QuadTo2m = generator.V2QuadTo2m() + V2mTo2Quad = generator.V2mTo2Quad() + + json_io.write_tensor(quadpoints, "quadpoints", filename) + json_io.write_tensor( + quadweights.reshape(-1, 1), "quadweights", filename + ) + json_io.write_tensor(resample, "resample", filename) + json_io.write_tensor(V2QuadTo2m, "V2QuadTo2m", filename) + json_io.write_tensor(V2mTo2Quad, "V2mTo2Quad", filename) + for a in range(0, 4): + for b in range(0, 4): + V3mTo2n = generator.V3mTo2n(a, b) + V3mTo2nTWDivM = generator.V3mTo2nTWDivM( + a, b, matorder=materialorder + ) + json_io.write_tensor(V3mTo2n, f"V3mTo2n({a},{b})", filename) + json_io.write_tensor( + V3mTo2nTWDivM, f"V3mTo2nTWDivM({a},{b})", filename + ) + elif sys.argv[1] == "hom": + materialorder = None + if len(sys.argv) > 2: + materialorder = int(sys.argv[2]) + for mode in ("nb", "ip"): + filename = f"hom-{mode}-h{materialorder}.json" + generator = PlasticityGenerator(materialorder) + points = generator.nodes(mode) + interpolator = generator.generate_Vandermonde_inv(mode) + deinterpolator = generator.generate_Vandermonde(mode) + + json_io.write_tensor(points, "hompoints", filename) + json_io.write_tensor(interpolator, "homproject", filename) + json_io.write_tensor(deinterpolator, "homiproject", filename) + for order in range(2, 9): + for quadname, quadrule in zip( + ["stroud", "dunavant", "witherden_vincent"], + [ + quad_points.stroud(order + 1), + quad_points.dunavant(order + 1), + quad_points.witherden_vincent(order + 1), + ], + ): + filename = f"homdr-{mode}-{quadname}-{order}-h{materialorder}.json" + generator = dr_generator(materialorder, quadrule) + for a in range(0, 4): + for b in range(0, 4): + V3mTo2n = generator.V3mTo2n(a, b) + json_io.write_tensor( + V3mTo2n, f"homV3mTo2n({a},{b})", filename + ) + + +if __name__ == "__main__": + main() diff --git a/src/seissol_matrices/nodal/nodes_1.txt b/src/seissol_matrices/nodal/nodes_1.txt new file mode 100644 index 0000000..baf2610 --- /dev/null +++ b/src/seissol_matrices/nodal/nodes_1.txt @@ -0,0 +1 @@ +0.33333333333333333333 0.33333333333333333333 diff --git a/src/seissol_matrices/plasticity/nb_1.txt b/src/seissol_matrices/plasticity/nb_1.txt new file mode 100644 index 0000000..abca404 --- /dev/null +++ b/src/seissol_matrices/plasticity/nb_1.txt @@ -0,0 +1 @@ +0.333333333333333 0.333333333333333 0.333333333333333 diff --git a/src/seissol_matrices/vtk_points.py b/src/seissol_matrices/vtk_points.py new file mode 100644 index 0000000..256dabd --- /dev/null +++ b/src/seissol_matrices/vtk_points.py @@ -0,0 +1,109 @@ +import numpy as np + + +def vtk_lagrange_2d(d, start=0, divd=None): + if d < 0: + return np.empty((0, 2)) + if d == 0: + return np.array([[1 / 3, 1 / 3]]) + else: + if divd is None: + divd = d + s = 1 / divd + base = 0 + obase = 0 + allrings = [] + for r in range(start, d // 3 + 1): + if d == 3 * r: + allrings += [[s * r, s * r]] + else: + ringcorner = [ + [s * r, s * r], + [s * (d - 2 * r), s * r], + [s * r, s * (d - 2 * r)], + ] + ringborder = [[s * (r + i), s * r] for i in range(1, d - 3 * r)] + ringborder += [ + [s * (d - 2 * r - i), s * (r + i)] for i in range(1, d - 3 * r) + ] + ringborder += [ + [s * r, s * (d - 2 * r - i)] for i in range(1, d - 3 * r) + ] + allrings += ringcorner + ringborder + return np.array(allrings) + + +def vtk_lagrange_3d(d): + if d == 0: + return np.array([[1 / 4, 1 / 4, 1 / 4]]) + else: + s = 1 / d + base = 0 + obase = 0 + allrings = [] + for r in range(d // 4 + 1): + if d == 4 * r: + allrings += [[s * r, s * r, s * r]] + else: + ringcorner = [ + [s * r, s * r, s * r], + [s * (d - 3 * r), s * r, s * r], + [s * r, s * (d - 3 * r), s * r], + [s * r, s * r, s * (d - 3 * r)], + ] + ringborder = [[s * (r + i), s * r, s * r] for i in range(1, d - 4 * r)] + ringborder += [ + [s * (d - 3 * r - i), s * (r + i), s * r] + for i in range(1, d - 4 * r) + ] + ringborder += [ + [s * r, s * (d - 3 * r - i), s * r] for i in range(1, d - 4 * r) + ] + ringborder += [[s * r, s * r, s * (r + i)] for i in range(1, d - 4 * r)] + ringborder += [ + [s * (d - 3 * r - i), s * r, s * (r + i)] + for i in range(1, d - 4 * r) + ] + ringborder += [ + [s * r, s * (d - 3 * r - i), s * (r + i)] + for i in range(1, d - 4 * r) + ] + subface = vtk_lagrange_2d(d - r, r + 1, d) + ringfaces = [(i, s * r, j) for i, j in subface] + ringfaces += [(j, (1 - s * r) - i - j, i) for i, j in subface] + ringfaces += [(s * r, j, i) for i, j in subface] + ringfaces += [(j, i, s * r) for i, j in subface] + allrings += ringcorner + ringborder + ringfaces + return np.array(allrings) + + +def vtk_lagrange_2d_from_vtk(d): + import vtk + + tet = vtk.vtkLagrangeTriangle() + points = ((d + 1) * (d + 2)) // 2 + tet.GetPointIds().SetNumberOfIds(points) + tet.GetPoints().SetNumberOfPoints(points) + tet.Initialize() + pointlist = [] + for i in range(points): + barycentric = [0, 0, 0] + tet.ToBarycentricIndex(i, barycentric) + pointlist += [[b / d for b in barycentric[:2]]] + return np.array(pointlist) + + +def vtk_lagrange_3d_from_vtk(d): + import vtk + + tet = vtk.vtkLagrangeTetra() + points = ((d + 1) * (d + 2) * (d + 3)) // 6 + tet.GetPointIds().SetNumberOfIds(points) + tet.GetPoints().SetNumberOfPoints(points) + tet.Initialize() + pointlist = [] + for i in range(points): + barycentric = [0, 0, 0, 0] + tet.ToBarycentricIndex(i, barycentric) + pointlist += [[b / d for b in barycentric[:3]]] + return np.array(pointlist)