NVIDIA Modulus Sym v1.0.0
Sym v1.0.0

# Source code for modulus.sym.utils.vpinn.integral

# Copyright (c) 2023, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

""" Helper functions and classes for integration
"""

import torch
import numpy as np

def tensor_int(w, v, u=None):
# u is a N*1 tensor
# v is a N*M tensor
# w is a N*1 tensor, quadrature (cubature) weights
# N is the number of points
# M is the number of (test) functions
# return: 1*M tensor: integrals of u*v[i] if u is not None
# return: 1*M tensor: integrals of v[i] if u is None
if u is None:
else:

# class template for the quadrature
def __init__(self, scheme, trans, jac):
# scheme is the quadpy scheme
# trans is the transform from reference domain to target domain
# jac is the jacobian of trans, SHOULD BE 1D NUMPY ARRAY!
# points_ref and weights_ref are on the reference domain
self.scheme = scheme
self.trans = trans
self.jac = jac

self.points_ref = scheme.points
self.weights_ref = scheme.weights
self.make_numpy()
# self.make_tensor()
self.N_points = self.points_numpy.shape[0]

def make_numpy(self):
# points_numpy and weights_numpy are N*d numpy arrays, where N is the # of points and d is the dimension
# The approximated integral value is given by np.dot(f(p[:,0],p[:,1],p[:,2]),w)
self.points_numpy = self.trans(self.points_ref)
self.weights_numpy = (
self.weights_ref * self.jac
)  # check here, should be 1D*1D numpy array, or 1D*constant

def make_tensor(self):
# points_tensor and weights_tensor are N*d tf tensors, where N is the # of points and d is the dimension
self.points_tensor = torch.tensor(self.points_numpy, dtype=torch.float32)
self.weights_tensor = torch.tensor(
self.weights_numpy.reshape((-1, 1)), dtype=torch.float32
)

def __init__(self, points_numpy, weights_numpy):
self.points_numpy = points_numpy
self.weights_numpy = weights_numpy
# self.make_tensor()

def make_tensor(self):
# points_tensor and weights_tensor are N*d tf tensors, where N is the # of points and d is the dimension
self.points_tensor = torch.tensor(self.points_numpy, dtype=torch.float32)
self.weights_tensor = torch.tensor(
self.weights_numpy.reshape((-1, 1)), dtype=torch.float32
)

points_tmp = []
weights_tmp = []
for para in paras:

# 1D classes. (Quad_Line can be used in nD)
def __init__(self, p0, p1, n, scheme_fcn=qd.c1.gauss_legendre):
self.p0 = np.reshape(np.array(p0), (1, -1))
self.p1 = np.reshape(np.array(p1), (1, -1))
super().__init__(
scheme=scheme_fcn(n),
trans=lambda t: 0.5 * (self.p0 + self.p1)
+ 0.5 * (self.p1 - self.p0) * np.reshape(t, (-1, 1)),
jac=np.linalg.norm(self.p1 - self.p0) / 2,
)

# 2D curves
def __init__(self, r, c, n, scheme_fcn=qd.u2.get_good_scheme):
self.r = np.array(r)
self.c = np.array(c)

def my_trans(x):
rr = np.multiply.outer(self.r, x)
rr = np.swapaxes(rr, 0, -2)
return rr + self.c

super().__init__(scheme=scheme_fcn(n), trans=my_trans, jac=2 * np.pi * self.r)

# 2D domains
def __init__(self, v, n, scheme_fcn=qd.t2.get_good_scheme):

self.v = np.array(v)  # 3x2 numpy array
if self.v.shape != (3, 2):
self.v = self.v.T
assert self.v.shape == (3, 2), "Vertices must be a 3 by 2 list or numpy array!"
self.vol = get_vol(self.v)
super().__init__(
scheme=scheme_fcn(n), trans=lambda x: x.T @ self.v, jac=self.vol
)

def __init__(self, r, c, n, scheme_fcn=qd.s2.get_good_scheme):
self.r = np.array(r)
self.c = np.array(c)

def my_trans(x):
rr = np.multiply.outer(self.r, x.T)
rr = np.swapaxes(rr, 0, -2)
return rr + self.c

super().__init__(scheme=scheme_fcn(n), trans=my_trans, jac=np.pi * self.r**2)

"""
The points are specified in an array of shape (2, 2, ...) such that arr[0][0] is the lower left corner, arr[1][1] the upper right, and set region_type=False.
If your c2 has its sides aligned with the coordinate axes, you can use v=[[x0, x1], [y0, y1]], and set region_type=True (default).
"""

def __init__(self, v, n, region_type=True, scheme_fcn=qd.c2.get_good_scheme):

if region_type:

self.v = rectangle_points(*v)
else:
self.v = v
super().__init__(
scheme=scheme_fcn(n),
trans=lambda x: transform(x, self.v),
jac=lambda x: np.abs(get_detJ(x, self.v))
* 2 ** np.prod(len(self.v.shape) - 1),
)

def make_numpy(self):
self.points_numpy = self.trans(self.points_ref)
self.weights_numpy = self.weights_ref * self.jac(
self.points_ref
)  # check here, should be 1D*1D numpy array, or 1D*constant# 3D surfaces
def __init__(self, r, c, n, scheme_fcn=qd.u3.get_good_scheme):
self.r = np.array(r)
self.c = np.array(c)
super().__init__(
scheme=scheme_fcn(n),
trans=lambda x: x.T * self.r + self.c,
jac=4 * np.pi * self.r**2,
)

# 3D domain
def __init__(self, r, c, n, scheme_fcn=qd.s3.get_good_scheme):
assert (
n <= 7
), "The degree of the cubature is not more than 7. Otherwise use nD ball scheme!"
self.r = np.array(r)
self.c = np.array(c)

def my_trans(x):
rr = np.multiply.outer(self.r, x.T)
rr = np.swapaxes(rr, 0, -2)
return rr + self.c

super().__init__(
scheme=scheme_fcn(n), trans=my_trans, jac=4 / 3 * np.pi * self.r**3
)

def __init__(self, v, n, scheme_fcn=qd.t3.get_good_scheme):
assert (
n <= 14
), "The degree of the cubature is not more than 14. Otherwise use nD simplex scheme!"
self.v = np.array(v)
if self.v.shape != (4, 3):
self.v = self.v.T
assert self.v.shape == (4, 3), "Vertices must be a 4 by 3 list or numpy array!"

self.vol = get_vol(self.v)
super().__init__(
scheme=scheme_fcn(n), trans=lambda x: transform(x, self.v.T).T, jac=self.vol
)

def __init__(self, v, n, region_type=True, scheme_fcn=qd.c3.get_good_scheme):

assert (
n <= 11
), "The degree of the cubature is not more than 11. Otherwise use nD cube scheme!"
if region_type:

self.v = cube_points(*v)
else:
self.v = v
super().__init__(
scheme=scheme_fcn(n),
trans=lambda x: transform(x, self.v),
jac=lambda x: np.abs(get_detJ(x, self.v))
* 2 ** np.prod(len(self.v.shape) - 1),
)

def make_numpy(self):
self.points_numpy = self.trans(self.points_ref)
self.weights_numpy = self.weights_ref * self.jac(
self.points_ref
)  # check here, should be 1D*1D numpy array, or 1D*constant

def __init__(self, v, scheme_fcn=qd.p3.felippa_5):

self.v = v
super().__init__(
scheme=scheme_fcn(),
trans=lambda x: _transform(x.T, self.v).T,
jac=lambda x: np.abs(_get_det_J(self.v, x.T)),
)

def make_numpy(self):
self.points_numpy = self.trans(self.points_ref)
self.weights_numpy = self.weights_ref * self.jac(
self.points_ref
)  # check here, should be 1D*1D numpy array, or 1D*constant

def __init__(self, v, scheme_fcn=qd.w3.felippa_6):

self.v = np.array(v)
super().__init__(
scheme=scheme_fcn(),
trans=lambda x: _transform(x.T, self.v).T,
jac=lambda x: np.abs(_get_detJ(x.T, self.v)),
)

def make_numpy(self):
self.points_numpy = self.trans(self.points_ref)
self.weights_numpy = self.weights_ref * self.jac(
self.points_ref
)  # check here, should be 1D*1D numpy array, or 1D*constant

# nD manifold
def __init__(self, r, c, dim, scheme_fcn=qd.un.dobrodeev_1978):
import ndim

self.r = np.array(r)
self.c = np.array(c)
self.dim = dim

def my_trans(x):
rr = np.multiply.outer(self.r, x)
rr = np.swapaxes(rr, 0, -2)
return rr + self.c

self.vol = ndim.nsphere.volume(self.dim, r=self.r)
super().__init__(scheme=scheme_fcn(self.dim), trans=my_trans, jac=self.vol)

def __init__(self, v, dim, n, scheme_fcn=qd.tn.grundmann_moeller):

self.v = np.array(v)
self.dim = dim
self.vol = get_vol(self.v)
super().__init__(
scheme=scheme_fcn(self.dim, n),
trans=lambda x: transform(x, self.v.T).T,
jac=self.vol,
)

def __init__(self, r, c, dim, scheme_fcn=qd.sn.dobrodeev_1970):
import ndim

self.r = np.array(r)
self.c = np.array(c)
self.dim = dim
self.vol = ndim.nball.volume(self.dim, r=self.r, symbolic=False)

def my_trans(x):
rr = np.multiply.outer(self.r, x.T)
rr = np.swapaxes(rr, 0, -2)
return rr + self.c

super().__init__(scheme=scheme_fcn(self.dim), trans=my_trans, jac=self.vol)

def __init__(self, v, dim, region_type=True, scheme_fcn=qd.cn.stroud_cn_3_3):

self.dim = dim
if region_type:

self.v = ncube_points(*v)
else:
self.v = v
super().__init__(
scheme=scheme_fcn(self.dim),
trans=lambda x: transform(x, self.v),
jac=lambda x: 2 ** np.prod(len(self.v.shape) - 1)
* np.abs(get_detJ(x, self.v)),
)

def make_numpy(self):
self.points_numpy = self.trans(self.points_ref)
self.weights_numpy = self.weights_ref * self.jac(
self.points_ref
)  # check here, should be 1D*1D numpy array, or 1D*constant

# 2D cubature based on mesh
def domain_weights_and_points_2D(P, T, n=5, scheme=None):
# P is the point info
# T is the triangle info
# n is the cubature order, if applicable
T = T.astype(np.int)
Nt = T.shape[0]
if scheme is None:
scheme = qd.t2._lether.lether(n)
p_ref = scheme.points
w_ref = scheme.weights
xy_tmp = []
w_tmp = []
for i in range(1, Nt):
idp = T[i, :]
tri = np.vstack((P[idp[0], :], P[idp[1], :], P[idp[2], :]))
S = 0.5 * np.abs(np.linalg.det(np.hstack((tri, np.ones((3, 1))))))
xy_tmp.append(p_ref.T @ tri)
w_tmp.append(S * w_ref)
xy = np.vstack(xy_tmp)
w = np.hstack(w_tmp)
return w.astype(np.float32), xy.astype(np.float32)

# 3D cubature based on mesh
def domain_weights_and_points_3D(P, T, n=5, scheme=None):
# P is the point info
# T is the triangle info
# n is the cubature order, if applicable
T = T.astype(np.int)
Nt = T.shape[0]
if scheme is None:
scheme = qd.t3.get_good_scheme(n)
p_ref = scheme.points
w_ref = scheme.weights
xyz_tmp = []
w_tmp = []
for i in range(0, Nt):
idp = T[i, :]
tet = np.vstack((P[idp[0], :], P[idp[1], :], P[idp[2], :], P[idp[3], :]))
V = np.abs(np.linalg.det(np.hstack((tet, np.ones((4, 1)))))) / 6
xyz_tmp.append(p_ref.T @ tet)
w_tmp.append(V * w_ref)
xyz = np.vstack(xyz_tmp)
w = np.hstack(w_tmp)
return w.astype(np.float32), xyz.astype(np.float32)

# Householder reflector
def Householder_reflector(u0, v0):
# u and v are unit vectors
# Hu=v, Hv=u
u = u0.reshape((-1, 1)) / np.linalg.norm(u0)
v = v0.reshape((-1, 1)) / np.linalg.norm(v0)
return np.eye(3) + (u @ v.T + v @ u.T - u @ u.T - v @ v.T) / (1 - u.T @ v)


© Copyright 2023, NVIDIA Modulus Team. Last updated on Aug 8, 2023.