# (C) Copyright IBM Corp. 2019, 2020, 2021, 2022.
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import numpy as np
from scipy.special import jacobi
from itertools import product
import sys
[docs]class GaussLegendre:
def __init__(self, p_order=None):
self.p_order = p_order
self.alpha = 0
self.beta = 0
self.reference_interval = [-1, 1]
# Considering a homogeneous p-order
# Evaluate the weights and nodes of the
# element
if isinstance(p_order, int):
self.poly = jacobi(self.p_order, self.alpha, self.beta)
self.poly_der = self.poly.deriv(1)
self.poly_roots = sorted(self.poly.roots)
self.exec = self._exec_homogeneous
self.weights = [2 / ((1 - root ** 2) * self.poly_der(root) ** 2)
for root in self.poly_roots]
print('')
elif isinstance(p_order, tuple):
self.poly = tuple()
self.poly_der = tuple()
self.poly_roots = tuple()
self.weights = tuple()
for _p_order in p_order:
poly = jacobi(_p_order, self.alpha, self.beta)
poly_der = poly.deriv(1)
poly_roots = sorted(poly.roots)
self.exec = self._execute_adaptative
weights = [2/((1 - root**2)*poly_der(root)**2)
for root in poly_roots]
self.poly += (poly,)
self.poly_der += (poly_der,)
self.poly_roots += (poly_roots,)
self.weights += (weights,)
else:
pass
def _exec_homogeneous(self):
pass
def _execute_adaptative(self):
pass
[docs] def generate_domain(self, mesh=None):
nodes = mesh.internal_product(self.poly_roots)
n_dim = mesh.n_dim
weights = np.array(mesh.internal_product(self.weights)).prod(axis=1)[:, None]
nodes_list = list()
weights_list = list()
for key, el in mesh.elements.items():
sys.stdout.write("\rMapping from the reference to the real mesh element {}".format(key))
sys.stdout.flush()
nodes_mapped = mesh.map_to_element(nodes, self.reference_interval, el)
nodes_list.append(nodes_mapped)
weights_list.append(weights)
nodes_array = np.vstack(nodes_list)
weights_array = np.vstack(weights_list)
dim_arrays = np.split(nodes_array, n_dim, axis=1)
dim_arrays.append(weights_array)
return tuple(dim_arrays)
[docs] def generate_boundaries(self, mesh=None):
boundaries_list = dict()
for boundary in mesh.boundary_nodes_tags:
nodes_list = list()
weights_list = list()
for key, (el, tag) in mesh.boundary_elements[boundary].items():
nodes = mesh.internal_boundary_product(self.poly_roots)
weights = np.array(mesh.internal_boundary_product(self.weights)).prod(axis=1)[:, None]
sys.stdout.write("\rMapping from the reference to the real mesh element {} from {}".format(key, boundary))
sys.stdout.flush()
if isinstance(self.p_order, tuple):
nodes_mapped = mesh.map_to_boundary_element(nodes, self.reference_interval, el, tag)
nodes_list.append(nodes_mapped)
weights_list.append(weights)
else:
nodes_mapped = mesh.map_to_boundary_element(nodes, self.reference_interval, el)
nodes_list.append(nodes_mapped.T)
weights_list.append(weights)
nodes_array = np.vstack(nodes_list)
weights_array = np.vstack(weights_list)
boundaries_list[boundary] = (nodes_array, weights_array)
return boundaries_list