From 550e8a155d369537ba96a5c7d7d3a56023a86d5e Mon Sep 17 00:00:00 2001 From: Sebastian Wolf Date: Thu, 25 Jan 2024 17:17:34 +0100 Subject: [PATCH 1/2] Add nb plasticity matrices --- src/seissol_matrices/plasticity.py | 51 ++++++++++ src/seissol_matrices/plasticity/nb_2.txt | 4 + src/seissol_matrices/plasticity/nb_3.txt | 10 ++ src/seissol_matrices/plasticity/nb_4.txt | 20 ++++ src/seissol_matrices/plasticity/nb_5.txt | 35 +++++++ src/seissol_matrices/plasticity/nb_6.txt | 56 +++++++++++ src/seissol_matrices/plasticity/nb_7.txt | 84 ++++++++++++++++ src/seissol_matrices/plasticity/nb_8.txt | 120 +++++++++++++++++++++++ tests/test_plasticity.py | 72 ++++++++++++++ 9 files changed, 452 insertions(+) create mode 100644 src/seissol_matrices/plasticity.py create mode 100644 src/seissol_matrices/plasticity/nb_2.txt create mode 100644 src/seissol_matrices/plasticity/nb_3.txt create mode 100644 src/seissol_matrices/plasticity/nb_4.txt create mode 100644 src/seissol_matrices/plasticity/nb_5.txt create mode 100644 src/seissol_matrices/plasticity/nb_6.txt create mode 100644 src/seissol_matrices/plasticity/nb_7.txt create mode 100644 src/seissol_matrices/plasticity/nb_8.txt create mode 100755 tests/test_plasticity.py diff --git a/src/seissol_matrices/plasticity.py b/src/seissol_matrices/plasticity.py new file mode 100644 index 0000000..e225270 --- /dev/null +++ b/src/seissol_matrices/plasticity.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python3 + +import numpy as np +import os +from seissol_matrices import basis_functions, dg_matrices, json_io + + +def parse_nodes(filename): + with open(filename) as f: + lines = [[float(e) for e in l.rstrip().split()] for l in f] + return np.array(lines) + + +class PlasticityGenerator: + def __init__(self, order): + self.bf3_generator = basis_functions.BasisFunctionGenerator3D(order) + self.dg3_generator = dg_matrices.dg_generator(order, d=3) + self.order = order + + def generate_Vandermonde(self, mode): + assert mode in ["ip", "nb"] + nodes = parse_nodes( + f"{os.path.dirname(__file__)}/plasticity/{mode}_{self.order}.txt" + ) + m = self.bf3_generator.number_of_basis_functions() + assert nodes.shape[0] == m + + vandermonde = np.zeros((m, m)) + + for i in range(m): + for j in range(m): + n = nodes[i] + vandermonde[i, j] = self.bf3_generator.eval_basis([n[0], n[1], n[2]], j) + + return vandermonde + + def generate_Vandermonde_inv(self, mode): + vandermonde = self.generate_Vandermonde(mode) + vandermonde_inv = np.linalg.solve(vandermonde, np.eye(vandermonde.shape[0])) + return vandermonde_inv + + +if __name__ == "__main__": + for mode in ["nb"]: + for order in range(2, 8): + generator = PlasticityGenerator(order) + nb_vandermonde = generator.generate_Vandermonde(mode) + nb_vandermonde_inv = generator.generate_Vandermonde_inv(mode) + filename = f"output/plasticity_{mode}_matrices_{order}.json" + xml_io.write_matrix(nb_vandermonde, "v", filename) + xml_io.write_matrix(nb_vandermonde_inv, "vInv", filename) diff --git a/src/seissol_matrices/plasticity/nb_2.txt b/src/seissol_matrices/plasticity/nb_2.txt new file mode 100644 index 0000000..bce0458 --- /dev/null +++ b/src/seissol_matrices/plasticity/nb_2.txt @@ -0,0 +1,4 @@ +0.000000000000000 0.000000000000000 0.000000000000000 +1.000000000000000 0.000000000000000 0.000000000000000 +0.000000000000000 1.000000000000000 0.000000000000000 +0.000000000000000 0.000000000000000 1.000000000000000 diff --git a/src/seissol_matrices/plasticity/nb_3.txt b/src/seissol_matrices/plasticity/nb_3.txt new file mode 100644 index 0000000..2b4689b --- /dev/null +++ b/src/seissol_matrices/plasticity/nb_3.txt @@ -0,0 +1,10 @@ + 0.000000000000000 0.000000000000000 0.000000000000000 + 0.500000000000000 0.000000000000000 0.000000000000000 + 1.000000000000000 0.000000000000000 0.000000000000000 + 0.000000000000000 0.500000000000000 0.000000000000000 + 0.500000000000000 0.500000000000000 0.000000000000000 + 0.000000000000000 1.000000000000000 0.000000000000000 + 0.000000000000000 0.000000000000000 0.500000000000000 + 0.500000000000000 0.000000000000000 0.500000000000000 + 0.000000000000000 0.500000000000000 0.500000000000000 + 0.000000000000000 0.000000000000000 1.000000000000000 diff --git a/src/seissol_matrices/plasticity/nb_4.txt b/src/seissol_matrices/plasticity/nb_4.txt new file mode 100644 index 0000000..c36307d --- /dev/null +++ b/src/seissol_matrices/plasticity/nb_4.txt @@ -0,0 +1,20 @@ + 0.000000000000000 0.000000000000000 0.000000000000000 + 0.276393202250021 0.000000000000000 0.000000000000000 + 0.723606797749979 0.000000000000000 0.000000000000000 + 1.000000000000000 0.000000000000000 0.000000000000000 + -0.000000000000000 0.276393202250021 0.000000000000000 + 0.333333333333333 0.333333333333333 0.000000000000000 + 0.723606797749979 0.276393202250021 0.000000000000000 + 0.000000000000000 0.723606797749979 0.000000000000000 + 0.276393202250021 0.723606797749979 0.000000000000000 + 0.000000000000000 1.000000000000000 0.000000000000000 + 0.000000000000000 0.000000000000000 0.276393202250021 + 0.333333333333333 0.000000000000000 0.333333333333333 + 0.723606797749979 0.000000000000000 0.276393202250021 + 0.000000000000000 0.333333333333333 0.333333333333333 + 0.333333333333333 0.333333333333333 0.333333333333333 + 0.000000000000000 0.723606797749979 0.276393202250021 + 0.000000000000000 0.000000000000000 0.723606797749979 + 0.276393202250021 0.000000000000000 0.723606797749979 + 0.000000000000000 0.276393202250021 0.723606797749979 + 0.000000000000000 0.000000000000000 1.000000000000000 diff --git a/src/seissol_matrices/plasticity/nb_5.txt b/src/seissol_matrices/plasticity/nb_5.txt new file mode 100644 index 0000000..8bf5b29 --- /dev/null +++ b/src/seissol_matrices/plasticity/nb_5.txt @@ -0,0 +1,35 @@ + 0.000000000000000 0.000000000000000 0.000000000000000 + 0.172673164646011 0.000000000000000 0.000000000000000 + 0.500000000000000 0.000000000000000 0.000000000000000 + 0.827326835353989 0.000000000000000 0.000000000000000 + 1.000000000000000 0.000000000000000 0.000000000000000 + 0.000000000000000 0.172673164646011 0.000000000000000 + 0.224208213954503 0.224208213954503 0.000000000000000 + 0.551583572090994 0.224208213954503 0.000000000000000 + 0.827326835353989 0.172673164646011 0.000000000000000 + 0.000000000000000 0.500000000000000 0.000000000000000 + 0.224208213954503 0.551583572090994 0.000000000000000 + 0.500000000000000 0.500000000000000 0.000000000000000 + 0.000000000000000 0.827326835353989 0.000000000000000 + 0.172673164646011 0.827326835353989 0.000000000000000 + 0.000000000000000 1.000000000000000 0.000000000000000 + 0.000000000000000 0.000000000000000 0.172673164646011 + 0.224208213954503 0.000000000000000 0.224208213954503 + 0.551583572090994 0.000000000000000 0.224208213954503 + 0.827326835353989 0.000000000000000 0.172673164646011 + 0.000000000000000 0.224208213954503 0.224208213954503 + 0.250000000000000 0.250000000000000 0.250000000000000 + 0.551583572090994 0.224208213954503 0.224208213954503 + 0.000000000000000 0.551583572090994 0.224208213954503 + 0.224208213954503 0.551583572090994 0.224208213954503 + 0.000000000000000 0.827326835353989 0.172673164646011 + 0.000000000000000 0.000000000000000 0.500000000000000 + 0.224208213954503 0.000000000000000 0.551583572090994 + 0.500000000000000 0.000000000000000 0.500000000000000 + 0.000000000000000 0.224208213954503 0.551583572090994 + 0.224208213954503 0.224208213954503 0.551583572090994 + 0.000000000000000 0.500000000000000 0.500000000000000 + -0.000000000000000 0.000000000000000 0.827326835353988 + 0.172673164646011 0.000000000000000 0.827326835353988 + 0.000000000000000 0.172673164646011 0.827326835353989 + 0.000000000000000 0.000000000000000 1.000000000000000 diff --git a/src/seissol_matrices/plasticity/nb_6.txt b/src/seissol_matrices/plasticity/nb_6.txt new file mode 100644 index 0000000..0a58c53 --- /dev/null +++ b/src/seissol_matrices/plasticity/nb_6.txt @@ -0,0 +1,56 @@ + 0.000000000000000 0.000000000000000 0.000000000000000 + 0.117472338035268 0.000000000000000 0.000000000000000 + 0.357384241759677 0.000000000000000 0.000000000000000 + 0.642615758240323 0.000000000000000 0.000000000000000 + 0.882527661964732 0.000000000000000 0.000000000000000 + 1.000000000000000 0.000000000000000 0.000000000000000 + -0.000000000000000 0.117472338035268 0.000000000000000 + 0.155728267675030 0.155728267675030 0.000000000000000 + 0.417123903159896 0.165752193680209 0.000000000000000 + 0.688543464649939 0.155728267675030 0.000000000000000 + 0.882527661964732 0.117472338035268 0.000000000000000 + 0.000000000000000 0.357384241759677 0.000000000000000 + 0.165752193680209 0.417123903159896 0.000000000000000 + 0.417123903159896 0.417123903159896 0.000000000000000 + 0.642615758240323 0.357384241759677 0.000000000000000 + 0.000000000000000 0.642615758240322 0.000000000000000 + 0.155728267675030 0.688543464649939 0.000000000000000 + 0.357384241759677 0.642615758240322 0.000000000000000 + 0.000000000000000 0.882527661964732 0.000000000000000 + 0.117472338035268 0.882527661964732 0.000000000000000 + 0.000000000000000 1.000000000000000 0.000000000000000 + 0.000000000000000 0.000000000000000 0.117472338035268 + 0.155728267675030 0.000000000000000 0.155728267675030 + 0.417123903159896 0.000000000000000 0.165752193680209 + 0.688543464649939 0.000000000000000 0.155728267675030 + 0.882527661964732 0.000000000000000 0.117472338035268 + -0.000000000000000 0.155728267675030 0.155728267675030 + 0.188834092390309 0.188834092390309 0.188834092390309 + 0.433497722829073 0.188834092390309 0.188834092390309 + 0.688543464649940 0.155728267675030 0.155728267675030 + -0.000000000000000 0.417123903159896 0.165752193680209 + 0.188834092390309 0.433497722829073 0.188834092390309 + 0.417123903159896 0.417123903159896 0.165752193680209 + 0.000000000000000 0.688543464649939 0.155728267675030 + 0.155728267675030 0.688543464649939 0.155728267675030 + 0.000000000000000 0.882527661964732 0.117472338035268 + -0.000000000000000 0.000000000000000 0.357384241759677 + 0.165752193680209 0.000000000000000 0.417123903159896 + 0.417123903159896 0.000000000000000 0.417123903159896 + 0.642615758240323 0.000000000000000 0.357384241759677 + 0.000000000000000 0.165752193680209 0.417123903159896 + 0.188834092390309 0.188834092390309 0.433497722829073 + 0.417123903159896 0.165752193680209 0.417123903159896 + 0.000000000000000 0.417123903159896 0.417123903159896 + 0.165752193680209 0.417123903159896 0.417123903159896 + 0.000000000000000 0.642615758240323 0.357384241759677 + 0.000000000000000 0.000000000000000 0.642615758240323 + 0.155728267675030 0.000000000000000 0.688543464649940 + 0.357384241759677 0.000000000000000 0.642615758240323 + 0.000000000000000 0.155728267675030 0.688543464649940 + 0.155728267675030 0.155728267675030 0.688543464649940 + 0.000000000000000 0.357384241759677 0.642615758240322 + 0.000000000000000 0.000000000000000 0.882527661964732 + 0.117472338035268 0.000000000000000 0.882527661964732 + 0.000000000000000 0.117472338035268 0.882527661964732 + 0.000000000000000 0.000000000000000 1.000000000000000 diff --git a/src/seissol_matrices/plasticity/nb_7.txt b/src/seissol_matrices/plasticity/nb_7.txt new file mode 100644 index 0000000..c0408cc --- /dev/null +++ b/src/seissol_matrices/plasticity/nb_7.txt @@ -0,0 +1,84 @@ + 0.000000000000000 0.000000000000000 0.000000000000000 + 0.084888051860716 0.000000000000000 0.000000000000000 + 0.265575603264643 0.000000000000000 0.000000000000000 + 0.500000000000000 0.000000000000000 0.000000000000000 + 0.734424396735357 0.000000000000000 0.000000000000000 + 0.915111948139283 0.000000000000000 0.000000000000000 + 1.000000000000000 0.000000000000000 0.000000000000000 + 0.000000000000000 0.084888051860716 0.000000000000000 + 0.113188538898858 0.113188538898858 0.000000000000000 + 0.319716439082216 0.120634457015483 0.000000000000000 + 0.559649103902302 0.120634457015483 0.000000000000000 + 0.773622922202284 0.113188538898858 0.000000000000000 + 0.915111948139284 0.084888051860716 0.000000000000000 + -0.000000000000000 0.265575603264643 0.000000000000000 + 0.120634457015483 0.319716439082216 0.000000000000000 + 0.333333333333333 0.333333333333333 0.000000000000000 + 0.559649103902302 0.319716439082216 0.000000000000000 + 0.734424396735357 0.265575603264643 0.000000000000000 + 0.000000000000000 0.500000000000000 0.000000000000000 + 0.120634457015483 0.559649103902302 0.000000000000000 + 0.319716439082216 0.559649103902302 0.000000000000000 + 0.500000000000000 0.500000000000000 0.000000000000000 + 0.000000000000000 0.734424396735357 0.000000000000000 + 0.113188538898858 0.773622922202284 0.000000000000000 + 0.265575603264643 0.734424396735357 0.000000000000000 + 0.000000000000000 0.915111948139283 0.000000000000000 + 0.084888051860717 0.915111948139283 0.000000000000000 + 0.000000000000000 1.000000000000000 0.000000000000000 + -0.000000000000000 0.000000000000000 0.084888051860717 + 0.113188538898858 0.000000000000000 0.113188538898858 + 0.319716439082216 0.000000000000000 0.120634457015483 + 0.559649103902302 0.000000000000000 0.120634457015483 + 0.773622922202284 0.000000000000000 0.113188538898858 + 0.915111948139283 0.000000000000000 0.084888051860717 + -0.000000000000000 0.113188538898858 0.113188538898858 + 0.144598596516112 0.144598596516112 0.144598596516112 + 0.347086254780990 0.152913745219009 0.152913745219009 + 0.566204210451665 0.144598596516112 0.144598596516112 + 0.773622922202284 0.113188538898858 0.113188538898858 + 0.000000000000000 0.319716439082216 0.120634457015483 + 0.152913745219009 0.347086254780991 0.152913745219009 + 0.347086254780991 0.347086254780991 0.152913745219009 + 0.559649103902302 0.319716439082216 0.120634457015483 + 0.000000000000000 0.559649103902302 0.120634457015483 + 0.144598596516112 0.566204210451665 0.144598596516112 + 0.319716439082216 0.559649103902302 0.120634457015483 + 0.000000000000000 0.773622922202284 0.113188538898858 + 0.113188538898858 0.773622922202284 0.113188538898858 + 0.000000000000000 0.915111948139284 0.084888051860716 + 0.000000000000000 0.000000000000000 0.265575603264643 + 0.120634457015483 0.000000000000000 0.319716439082216 + 0.333333333333333 0.000000000000000 0.333333333333333 + 0.559649103902302 0.000000000000000 0.319716439082216 + 0.734424396735357 0.000000000000000 0.265575603264643 + 0.000000000000000 0.120634457015483 0.319716439082216 + 0.152913745219009 0.152913745219009 0.347086254780991 + 0.347086254780991 0.152913745219009 0.347086254780991 + 0.559649103902302 0.120634457015483 0.319716439082216 + 0.000000000000000 0.333333333333333 0.333333333333333 + 0.152913745219009 0.347086254780991 0.347086254780991 + 0.333333333333333 0.333333333333333 0.333333333333333 + 0.000000000000000 0.559649103902302 0.319716439082216 + 0.120634457015483 0.559649103902302 0.319716439082216 + 0.000000000000000 0.734424396735357 0.265575603264643 + 0.000000000000000 0.000000000000000 0.500000000000000 + 0.120634457015483 0.000000000000000 0.559649103902302 + 0.319716439082216 0.000000000000000 0.559649103902302 + 0.500000000000000 0.000000000000000 0.500000000000000 + -0.000000000000000 0.120634457015483 0.559649103902302 + 0.144598596516112 0.144598596516112 0.566204210451665 + 0.319716439082216 0.120634457015483 0.559649103902302 + 0.000000000000000 0.319716439082216 0.559649103902302 + 0.120634457015483 0.319716439082216 0.559649103902302 + 0.000000000000000 0.500000000000000 0.500000000000000 + 0.000000000000000 0.000000000000000 0.734424396735357 + 0.113188538898858 0.000000000000000 0.773622922202284 + 0.265575603264643 0.000000000000000 0.734424396735357 + 0.000000000000000 0.113188538898858 0.773622922202284 + 0.113188538898858 0.113188538898858 0.773622922202284 + 0.000000000000000 0.265575603264643 0.734424396735357 + -0.000000000000000 0.000000000000000 0.915111948139284 + 0.084888051860716 0.000000000000000 0.915111948139284 + 0.000000000000000 0.084888051860717 0.915111948139283 + 0.000000000000000 0.000000000000000 1.000000000000000 diff --git a/src/seissol_matrices/plasticity/nb_8.txt b/src/seissol_matrices/plasticity/nb_8.txt new file mode 100644 index 0000000..857c5d3 --- /dev/null +++ b/src/seissol_matrices/plasticity/nb_8.txt @@ -0,0 +1,120 @@ + 0.000000000000000 0.000000000000000 0.000000000000000 + 0.064129925745196 0.000000000000000 0.000000000000000 + 0.204149909283429 0.000000000000000 0.000000000000000 + 0.395350391048761 0.000000000000000 0.000000000000000 + 0.604649608951239 0.000000000000000 0.000000000000000 + 0.795850090716571 0.000000000000000 0.000000000000000 + 0.935870074254803 0.000000000000000 0.000000000000000 + 1.000000000000000 0.000000000000000 0.000000000000000 + -0.000000000000000 0.064129925745197 0.000000000000000 + 0.087358091648967 0.087358091648967 0.000000000000000 + 0.248894314819442 0.096650347272693 0.000000000000000 + 0.450520851832608 0.098958296334783 0.000000000000000 + 0.654455337907865 0.096650347272693 0.000000000000000 + 0.825283816702066 0.087358091648967 0.000000000000000 + 0.935870074254803 0.064129925745197 0.000000000000000 + 0.000000000000000 0.204149909283429 0.000000000000000 + 0.096650347272693 0.248894314819442 0.000000000000000 + 0.266664277434014 0.266664277434014 0.000000000000000 + 0.466671445131973 0.266664277434014 0.000000000000000 + 0.654455337907865 0.248894314819442 0.000000000000000 + 0.795850090716571 0.204149909283429 0.000000000000000 + 0.000000000000000 0.395350391048761 0.000000000000000 + 0.098958296334783 0.450520851832608 0.000000000000000 + 0.266664277434014 0.466671445131973 0.000000000000000 + 0.450520851832608 0.450520851832608 0.000000000000000 + 0.604649608951239 0.395350391048761 0.000000000000000 + 0.000000000000000 0.604649608951239 0.000000000000000 + 0.096650347272693 0.654455337907865 0.000000000000000 + 0.248894314819442 0.654455337907865 0.000000000000000 + 0.395350391048761 0.604649608951239 0.000000000000000 + 0.000000000000000 0.795850090716571 0.000000000000000 + 0.087358091648967 0.825283816702066 0.000000000000000 + 0.204149909283429 0.795850090716571 0.000000000000000 + 0.000000000000000 0.935870074254803 0.000000000000000 + 0.064129925745197 0.935870074254803 0.000000000000000 + 0.000000000000000 1.000000000000000 0.000000000000000 + 0.000000000000000 -0.000000000000000 0.064129925745197 + 0.087358091648967 0.000000000000000 0.087358091648967 + 0.248894314819442 0.000000000000000 0.096650347272693 + 0.450520851832608 0.000000000000000 0.098958296334783 + 0.654455337907865 0.000000000000000 0.096650347272693 + 0.825283816702066 0.000000000000000 0.087358091648967 + 0.935870074254803 0.000000000000000 0.064129925745197 + 0.000000000000000 0.087358091648967 0.087358091648967 + 0.115151061954242 0.115151061954242 0.115151061954242 + 0.277770654126193 0.126022787099008 0.126022787099008 + 0.470183771675792 0.126022787099008 0.126022787099008 + 0.654546814137274 0.115151061954242 0.115151061954242 + 0.825283816702066 0.087358091648967 0.087358091648967 + 0.000000000000000 0.248894314819442 0.096650347272693 + 0.126022787099008 0.277770654126193 0.126022787099008 + 0.288950992577170 0.288950992577170 0.133147022268490 + 0.470183771675792 0.277770654126193 0.126022787099008 + 0.654455337907865 0.248894314819442 0.096650347272693 + 0.000000000000000 0.450520851832608 0.098958296334783 + 0.126022787099008 0.470183771675792 0.126022787099008 + 0.277770654126193 0.470183771675792 0.126022787099008 + 0.450520851832608 0.450520851832608 0.098958296334783 + 0.000000000000000 0.654455337907865 0.096650347272693 + 0.115151061954242 0.654546814137274 0.115151061954242 + 0.248894314819442 0.654455337907865 0.096650347272693 + 0.000000000000000 0.825283816702066 0.087358091648967 + 0.087358091648967 0.825283816702066 0.087358091648967 + 0.000000000000000 0.935870074254803 0.064129925745197 + 0.000000000000000 0.000000000000000 0.204149909283429 + 0.096650347272693 0.000000000000000 0.248894314819442 + 0.266664277434013 0.000000000000000 0.266664277434014 + 0.466671445131973 0.000000000000000 0.266664277434014 + 0.654455337907865 0.000000000000000 0.248894314819442 + 0.795850090716571 0.000000000000000 0.204149909283429 + 0.000000000000000 0.096650347272693 0.248894314819442 + 0.126022787099008 0.126022787099008 0.277770654126193 + 0.288950992577170 0.133147022268490 0.288950992577170 + 0.470183771675792 0.126022787099008 0.277770654126193 + 0.654455337907865 0.096650347272693 0.248894314819442 + 0.000000000000000 0.266664277434014 0.266664277434014 + 0.133147022268490 0.288950992577170 0.288950992577170 + 0.288950992577170 0.288950992577170 0.288950992577170 + 0.466671445131973 0.266664277434014 0.266664277434014 + 0.000000000000000 0.466671445131973 0.266664277434014 + 0.126022787099008 0.470183771675792 0.277770654126193 + 0.266664277434014 0.466671445131973 0.266664277434014 + 0.000000000000000 0.654455337907865 0.248894314819442 + 0.096650347272693 0.654455337907865 0.248894314819442 + 0.000000000000000 0.795850090716571 0.204149909283429 + 0.000000000000000 0.000000000000000 0.395350391048761 + 0.098958296334783 0.000000000000000 0.450520851832608 + 0.266664277434014 0.000000000000000 0.466671445131973 + 0.450520851832608 0.000000000000000 0.450520851832608 + 0.604649608951239 0.000000000000000 0.395350391048761 + 0.000000000000000 0.098958296334783 0.450520851832608 + 0.126022787099008 0.126022787099008 0.470183771675792 + 0.277770654126193 0.126022787099008 0.470183771675792 + 0.450520851832608 0.098958296334783 0.450520851832608 + 0.000000000000000 0.266664277434014 0.466671445131973 + 0.126022787099008 0.277770654126193 0.470183771675792 + 0.266664277434014 0.266664277434014 0.466671445131973 + 0.000000000000000 0.450520851832608 0.450520851832608 + 0.098958296334783 0.450520851832608 0.450520851832608 + 0.000000000000000 0.604649608951239 0.395350391048761 + 0.000000000000000 0.000000000000000 0.604649608951239 + 0.096650347272693 0.000000000000000 0.654455337907865 + 0.248894314819442 0.000000000000000 0.654455337907865 + 0.395350391048761 0.000000000000000 0.604649608951239 + 0.000000000000000 0.096650347272693 0.654455337907865 + 0.115151061954242 0.115151061954242 0.654546814137274 + 0.248894314819442 0.096650347272693 0.654455337907865 + -0.000000000000000 0.248894314819442 0.654455337907865 + 0.096650347272693 0.248894314819442 0.654455337907865 + 0.000000000000000 0.395350391048761 0.604649608951239 + 0.000000000000000 0.000000000000000 0.795850090716571 + 0.087358091648967 0.000000000000000 0.825283816702066 + 0.204149909283429 0.000000000000000 0.795850090716571 + 0.000000000000000 0.087358091648967 0.825283816702066 + 0.087358091648967 0.087358091648967 0.825283816702066 + -0.000000000000000 0.204149909283429 0.795850090716571 + -0.000000000000000 0.000000000000000 0.935870074254803 + 0.064129925745197 0.000000000000000 0.935870074254803 + 0.000000000000000 0.064129925745197 0.935870074254803 + 0.000000000000000 0.000000000000000 1.000000000000000 diff --git a/tests/test_plasticity.py b/tests/test_plasticity.py new file mode 100755 index 0000000..ce4df54 --- /dev/null +++ b/tests/test_plasticity.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 + +import unittest +import numpy as np +import os +import requests + +from seissol_matrices import plasticity +from seissol_matrices import xml_io +from seissol_matrices import quad_points +from tests import helper + + +class abstract_tester(object): + def test_nb(self): + vandermonde = self.generator.generate_Vandermonde("nb") + vandermonde_from_file = xml_io.read_matrix("v", self.filename) + self.compare_matrices(vandermonde, vandermonde_from_file) + + vandermonde_inv = self.generator.generate_Vandermonde_inv("nb") + vandermonde_inv_from_file = xml_io.read_matrix("vInv", self.filename) + self.compare_matrices(vandermonde_inv, vandermonde_inv_from_file) + + +def setUpClassFromOrder(cls, order): + cls.filename = f"plasticity_nb_{order}.xml" + # point to some old commit + url = f"https://raw.githubusercontent.com/SeisSol/SeisSol/c99ca8f814b7184eb435272f9d4f63b03b8b6cf4/generated_code/matrices/plasticity_nb_matrices_{order}.xml" + r = requests.get(url, allow_redirects=True) + open(cls.filename, "wb").write(r.content) + cls.order = order + cls.generator = plasticity.PlasticityGenerator(cls.order) + + +class test_plasticity_2(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, 2) + + +class test_plasticity_3(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, 3) + + +class test_plasticity_4(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, 4) + + +class test_plasticity_5(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, 5) + + +class test_plasticity_6(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, 6) + + +class test_plasticity_7(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, 7) + + +if __name__ == "__main__": + unittest.main() From 6f9c7aa5096e2c1ac4d35d0b9bdcd2eb39837afe Mon Sep 17 00:00:00 2001 From: Sebastian Wolf Date: Fri, 26 Jan 2024 12:02:30 +0100 Subject: [PATCH 2/2] Add IP matrices --- src/seissol_matrices/plasticity.py | 51 ------------ src/seissol_matrices/plasticity/stroud.py | 41 ++++++++++ src/seissol_matrices/plasticity_matrices.py | 66 ++++++++++++++++ tests/test_plasticity.py | 88 ++++++++++++++++----- 4 files changed, 174 insertions(+), 72 deletions(-) delete mode 100644 src/seissol_matrices/plasticity.py create mode 100644 src/seissol_matrices/plasticity/stroud.py create mode 100644 src/seissol_matrices/plasticity_matrices.py diff --git a/src/seissol_matrices/plasticity.py b/src/seissol_matrices/plasticity.py deleted file mode 100644 index e225270..0000000 --- a/src/seissol_matrices/plasticity.py +++ /dev/null @@ -1,51 +0,0 @@ -#!/usr/bin/env python3 - -import numpy as np -import os -from seissol_matrices import basis_functions, dg_matrices, json_io - - -def parse_nodes(filename): - with open(filename) as f: - lines = [[float(e) for e in l.rstrip().split()] for l in f] - return np.array(lines) - - -class PlasticityGenerator: - def __init__(self, order): - self.bf3_generator = basis_functions.BasisFunctionGenerator3D(order) - self.dg3_generator = dg_matrices.dg_generator(order, d=3) - self.order = order - - def generate_Vandermonde(self, mode): - assert mode in ["ip", "nb"] - nodes = parse_nodes( - f"{os.path.dirname(__file__)}/plasticity/{mode}_{self.order}.txt" - ) - m = self.bf3_generator.number_of_basis_functions() - assert nodes.shape[0] == m - - vandermonde = np.zeros((m, m)) - - for i in range(m): - for j in range(m): - n = nodes[i] - vandermonde[i, j] = self.bf3_generator.eval_basis([n[0], n[1], n[2]], j) - - return vandermonde - - def generate_Vandermonde_inv(self, mode): - vandermonde = self.generate_Vandermonde(mode) - vandermonde_inv = np.linalg.solve(vandermonde, np.eye(vandermonde.shape[0])) - return vandermonde_inv - - -if __name__ == "__main__": - for mode in ["nb"]: - for order in range(2, 8): - generator = PlasticityGenerator(order) - nb_vandermonde = generator.generate_Vandermonde(mode) - nb_vandermonde_inv = generator.generate_Vandermonde_inv(mode) - filename = f"output/plasticity_{mode}_matrices_{order}.json" - xml_io.write_matrix(nb_vandermonde, "v", filename) - xml_io.write_matrix(nb_vandermonde_inv, "vInv", filename) diff --git a/src/seissol_matrices/plasticity/stroud.py b/src/seissol_matrices/plasticity/stroud.py new file mode 100644 index 0000000..962b353 --- /dev/null +++ b/src/seissol_matrices/plasticity/stroud.py @@ -0,0 +1,41 @@ +import sys + +sys.path.append("../..") +from quad_rules import GaussJacobi +import numpy as np + +np.set_printoptions(precision=15, floatmode="fixed") + + +def stroud(n): + points_0, weights_0 = GaussJacobi.GaussJacobi(0, 0).compute_nodes_and_weights(n) + points_1, weights_1 = GaussJacobi.GaussJacobi(1, 0).compute_nodes_and_weights(n) + points_2, weights_2 = GaussJacobi.GaussJacobi(2, 0).compute_nodes_and_weights(n) + + x = np.zeros((n**3, 3)) + w = np.zeros(n**3) + for i, index in enumerate(np.ndindex(n, n, n)): + x[i, :] = np.array([points_2[index[0]], points_1[index[1]], points_0[index[2]]]) + w[i] = ( + weights_2[index[0]] * weights_1[index[1]] * weights_0[index[2]] * 0.5**6 + ) + x = 0.5 * x + 0.5 + + y = np.zeros((n**3, 3)) + y[:, 0] = x[:, 0] + y[:, 1] = x[:, 1] * (1 - y[:, 0]) + y[:, 2] = x[:, 2] * (1 - y[:, 0] - y[:, 1]) + + return y, w + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser() + parser.add_argument("n", type=int) + args = parser.parse_args() + + y, w = stroud(args.n) + for p in y: + print(f"{p[0]} {p[1]} {p[2]}") diff --git a/src/seissol_matrices/plasticity_matrices.py b/src/seissol_matrices/plasticity_matrices.py new file mode 100644 index 0000000..e3823ce --- /dev/null +++ b/src/seissol_matrices/plasticity_matrices.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python3 + +import numpy as np +import os +from seissol_matrices import basis_functions, dg_matrices, xml_io +from seissol_matrices.plasticity import stroud + + +def parse_nodes(filename): + with open(filename) as f: + lines = [[float(e) for e in l.rstrip().split()] for l in f] + return np.array(lines) + + +class PlasticityGenerator: + def __init__(self, order): + self.bf3_generator = basis_functions.BasisFunctionGenerator3D(order) + self.dg3_generator = dg_matrices.dg_generator(order, d=3) + self.order = order + + def generate_Vandermonde(self, mode): + assert mode in ["ip", "nb"] + if mode == "nb": + nodes = parse_nodes( + f"{os.path.dirname(__file__)}/plasticity/nb_{self.order}.txt" + ) + else: + nodes, _ = stroud.stroud(self.order + 1) + m = self.bf3_generator.number_of_basis_functions() + n = nodes.shape[0] + + vandermonde = np.zeros((n, m)) + + for i in range(n): + for j in range(m): + node = nodes[i] + vandermonde[i, j] = self.bf3_generator.eval_basis( + [node[0], node[1], node[2]], j + ) + + return vandermonde + + def generate_Vandermonde_inv(self, mode): + vandermonde = self.generate_Vandermonde(mode) + if mode == "nb": + vandermonde_inv = np.linalg.solve(vandermonde, np.eye(vandermonde.shape[0])) + else: + mass = self.dg3_generator.mass_matrix() + nodes, weights = stroud.stroud(self.order + 1) + vandermonde_inv = np.zeros(vandermonde.T.shape) + for l in range(vandermonde_inv.shape[0]): + for i in range(vandermonde_inv.shape[1]): + vandermonde_inv[l, i] = vandermonde[i, l] * weights[i] / mass[l, l] + + return vandermonde_inv + + +if __name__ == "__main__": + for mode in ["nb", "ip"]: + for order in range(2, 8): + generator = PlasticityGenerator(order) + vandermonde = generator.generate_Vandermonde(mode) + vandermonde_inv = generator.generate_Vandermonde_inv(mode) + filename = f"output/plasticity_{mode}_matrices_{order}.xml" + xml_io.write_matrix(vandermonde, "v", filename) + xml_io.write_matrix(vandermonde_inv, "vInv", filename) diff --git a/tests/test_plasticity.py b/tests/test_plasticity.py index ce4df54..1fe2ae5 100755 --- a/tests/test_plasticity.py +++ b/tests/test_plasticity.py @@ -5,67 +5,113 @@ import os import requests -from seissol_matrices import plasticity +from seissol_matrices import plasticity_matrices from seissol_matrices import xml_io from seissol_matrices import quad_points from tests import helper +def filename(mode, order): + return f"plasticity_{mode}_{order}.xml" + + class abstract_tester(object): - def test_nb(self): - vandermonde = self.generator.generate_Vandermonde("nb") + def test(self): + vandermonde = self.generator.generate_Vandermonde(self.mode) vandermonde_from_file = xml_io.read_matrix("v", self.filename) self.compare_matrices(vandermonde, vandermonde_from_file) - vandermonde_inv = self.generator.generate_Vandermonde_inv("nb") + vandermonde_inv = self.generator.generate_Vandermonde_inv(self.mode) vandermonde_inv_from_file = xml_io.read_matrix("vInv", self.filename) self.compare_matrices(vandermonde_inv, vandermonde_inv_from_file) -def setUpClassFromOrder(cls, order): - cls.filename = f"plasticity_nb_{order}.xml" +def downlaod_matrix(mode, order): # point to some old commit - url = f"https://raw.githubusercontent.com/SeisSol/SeisSol/c99ca8f814b7184eb435272f9d4f63b03b8b6cf4/generated_code/matrices/plasticity_nb_matrices_{order}.xml" + url = f"https://raw.githubusercontent.com/SeisSol/SeisSol/c99ca8f814b7184eb435272f9d4f63b03b8b6cf4/generated_code/matrices/plasticity_{mode}_matrices_{order}.xml" r = requests.get(url, allow_redirects=True) - open(cls.filename, "wb").write(r.content) + fn = filename(mode, order) + open(fn, "wb").write(r.content) + + +def setUpClassFromOrder(cls, mode, order): + downlaod_matrix(mode, order) + cls.filename = filename(mode, order) + cls.mode = mode cls.order = order - cls.generator = plasticity.PlasticityGenerator(cls.order) + cls.generator = plasticity_matrices.PlasticityGenerator(cls.order) + + +class test_nb_2(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, "nb", 2) + + +class test_nb_3(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, "nb", 3) + + +class test_nb_4(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, "nb", 4) + + +class test_nb_5(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, "nb", 5) + + +class test_nb_6(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, "nb", 6) + + +class test_nb_7(abstract_tester, helper.helper): + @classmethod + def setUpClass(cls): + setUpClassFromOrder(cls, "nb", 7) -class test_plasticity_2(abstract_tester, helper.helper): +class test_ip_2(abstract_tester, helper.helper): @classmethod def setUpClass(cls): - setUpClassFromOrder(cls, 2) + setUpClassFromOrder(cls, "ip", 2) -class test_plasticity_3(abstract_tester, helper.helper): +class test_ip_3(abstract_tester, helper.helper): @classmethod def setUpClass(cls): - setUpClassFromOrder(cls, 3) + setUpClassFromOrder(cls, "ip", 3) -class test_plasticity_4(abstract_tester, helper.helper): +class test_ip_4(abstract_tester, helper.helper): @classmethod def setUpClass(cls): - setUpClassFromOrder(cls, 4) + setUpClassFromOrder(cls, "ip", 4) -class test_plasticity_5(abstract_tester, helper.helper): +class test_ip_5(abstract_tester, helper.helper): @classmethod def setUpClass(cls): - setUpClassFromOrder(cls, 5) + setUpClassFromOrder(cls, "ip", 5) -class test_plasticity_6(abstract_tester, helper.helper): +class test_ip_6(abstract_tester, helper.helper): @classmethod def setUpClass(cls): - setUpClassFromOrder(cls, 6) + setUpClassFromOrder(cls, "ip", 6) -class test_plasticity_7(abstract_tester, helper.helper): +class test_ip_7(abstract_tester, helper.helper): @classmethod def setUpClass(cls): - setUpClassFromOrder(cls, 7) + setUpClassFromOrder(cls, "ip", 7) if __name__ == "__main__":