Unverified Commit a2a85975 authored by Felix Schindler's avatar Felix Schindler Committed by GitHub

Merge pull request #1004 from pymor/update-builtin-discretizer

Refactor builtin cg discretizer
parents ca093d87 f6f8eac0
Pipeline #65721 passed with stages
in 47 minutes and 23 seconds
......@@ -24,6 +24,31 @@ from pymor.operators.numpy import NumpyMatrixBasedOperator
from pymor.vectorarrays.numpy import NumpyVectorSpace
LagrangeShapeFunctions = {
line: {1: [lambda X: 1 - X[..., 0],
lambda X: X[..., 0]]},
square: {1: [lambda X: (1 - X[..., 0]) * (1 - X[..., 1]),
lambda X: (1 - X[..., 1]) * (X[..., 0]),
lambda X: (X[..., 0]) * (X[..., 1]),
lambda X: (X[..., 1]) * (1 - X[..., 0])]},
triangle: {1: [lambda X: 1 - X[..., 0] - X[..., 1],
lambda X: X[..., 0],
lambda X: X[..., 1]]},
}
LagrangeShapeFunctionsGrads = {
line: {1: np.array(([-1.],
[1., ]))},
square: {1: lambda X: np.array(([X[..., 1] - 1., X[..., 0] - 1.],
[1. - X[..., 1], - X[..., 0]],
[X[..., 1], X[..., 0]],
[-X[..., 1], 1. - X[..., 0]]))},
triangle: {1: np.array(([-1., -1.],
[1., 0.],
[0., 1.]))}
}
def CGVectorSpace(grid, id='STATE'):
return NumpyVectorSpace(grid.size(grid.dim), id)
......@@ -66,14 +91,8 @@ class L2ProductFunctionalP1(NumpyMatrixBasedOperator):
# evaluate the shape functions at the quadrature points on the reference
# element -> shape = (number of shape functions, number of quadrature points)
q, w = g.reference_element.quadrature(order=1)
if g.dim == 1:
SF = np.array((1 - q[..., 0], q[..., 0]))
elif g.dim == 2:
SF = np.array(((1 - np.sum(q, axis=-1)),
q[..., 0],
q[..., 1]))
else:
raise NotImplementedError
SF = LagrangeShapeFunctions[g.reference_element][1]
SF = np.array(tuple(f(q) for f in SF))
# integrate the products of the function with the shape functions on each element
# -> shape = (g.size(0), number of shape functions)
......@@ -219,13 +238,8 @@ class L2ProductFunctionalQ1(NumpyMatrixBasedOperator):
# evaluate the shape functions at the quadrature points on the reference
# element -> shape = (number of shape functions, number of quadrature points)
q, w = g.reference_element.quadrature(order=1)
if g.dim == 2:
SF = np.array(((1 - q[..., 0]) * (1 - q[..., 1]),
(1 - q[..., 1]) * (q[..., 0]),
(q[..., 0]) * (q[..., 1]),
(q[..., 1]) * (1 - q[..., 0])))
else:
raise NotImplementedError
SF = LagrangeShapeFunctions[g.reference_element][1]
SF = np.array(tuple(f(q) for f in SF))
# integrate the products of the function with the shape functions on each element
# -> shape = (g.size(0), number of shape functions)
......@@ -283,32 +297,21 @@ class L2ProductP1(NumpyMatrixBasedOperator):
g = self.grid
bi = self.boundary_info
# our shape functions
if g.dim == 2:
SF = [lambda X: 1 - X[..., 0] - X[..., 1],
lambda X: X[..., 0],
lambda X: X[..., 1]]
q, w = triangle.quadrature(order=2)
elif g.dim == 1:
SF = [lambda X: 1 - X[..., 0],
lambda X: X[..., 0]]
q, w = line.quadrature(order=2)
else:
raise NotImplementedError
# evaluate the shape functions on the quadrature points
SFQ = np.array(tuple(f(q) for f in SF))
q, w = g.reference_element.quadrature(order=2)
SF = LagrangeShapeFunctions[g.reference_element][1]
SF = np.array(tuple(f(q) for f in SF))
self.logger.info('Integrate the products of the shape functions on each element')
# -> shape = (g.size(0), number of shape functions ** 2)
if self.coefficient_function is not None:
C = self.coefficient_function(self.grid.centers(0), mu=mu)
SF_INTS = np.einsum('iq,jq,q,e,e->eij', SFQ, SFQ, w, g.integration_elements(0), C).ravel()
SF_INTS = np.einsum('iq,jq,q,e,e->eij', SF, SF, w, g.integration_elements(0), C).ravel()
del C
else:
SF_INTS = np.einsum('iq,jq,q,e->eij', SFQ, SFQ, w, g.integration_elements(0)).ravel()
SF_INTS = np.einsum('iq,jq,q,e->eij', SF, SF, w, g.integration_elements(0)).ravel()
del SFQ
del SF
self.logger.info('Determine global dofs ...')
SF_I0 = np.repeat(g.subentities(0, g.dim), g.dim + 1, axis=1).ravel()
......@@ -373,30 +376,21 @@ class L2ProductQ1(NumpyMatrixBasedOperator):
g = self.grid
bi = self.boundary_info
# our shape functions
if g.dim == 2:
SF = [lambda X: (1 - X[..., 0]) * (1 - X[..., 1]),
lambda X: (1 - X[..., 1]) * (X[..., 0]),
lambda X: (X[..., 0]) * (X[..., 1]),
lambda X: (1 - X[..., 0]) * (X[..., 1])]
else:
raise NotImplementedError
q, w = square.quadrature(order=2)
# evaluate the shape functions on the quadrature points
SFQ = np.array(tuple(f(q) for f in SF))
q, w = square.quadrature(order=2)
SF = LagrangeShapeFunctions[g.reference_element][1]
SF = np.array(tuple(f(q) for f in SF))
self.logger.info('Integrate the products of the shape functions on each element')
# -> shape = (g.size(0), number of shape functions ** 2)
if self.coefficient_function is not None:
C = self.coefficient_function(self.grid.centers(0), mu=mu)
SF_INTS = np.einsum('iq,jq,q,e,e->eij', SFQ, SFQ, w, g.integration_elements(0), C).ravel()
SF_INTS = np.einsum('iq,jq,q,e,e->eij', SF, SF, w, g.integration_elements(0), C).ravel()
del C
else:
SF_INTS = np.einsum('iq,jq,q,e->eij', SFQ, SFQ, w, g.integration_elements(0)).ravel()
SF_INTS = np.einsum('iq,jq,q,e->eij', SF, SF, w, g.integration_elements(0)).ravel()
del SFQ
del SF
self.logger.info('Determine global dofs ...')
SF_I0 = np.repeat(g.subentities(0, g.dim), 4, axis=1).ravel()
......@@ -474,17 +468,9 @@ class DiffusionOperatorP1(NumpyMatrixBasedOperator):
bi = self.boundary_info
# gradients of shape functions
if g.dim == 2:
SF_GRAD = np.array(([-1., -1.],
[1., 0.],
[0., 1.]))
elif g.dim == 1:
SF_GRAD = np.array(([-1.],
[1., ]))
else:
raise NotImplementedError
SF_GRAD = LagrangeShapeFunctionsGrads[g.reference_element][1]
self.logger.info('Calulate gradients of shape functions transformed by reference map ...')
self.logger.info('Calculate gradients of shape functions transformed by reference map ...')
SF_GRADS = np.einsum('eij,pj->epi', g.jacobian_inverse_transposed(0), SF_GRAD)
self.logger.info('Calculate all local scalar products between gradients ...')
......@@ -586,19 +572,14 @@ class DiffusionOperatorQ1(NumpyMatrixBasedOperator):
bi = self.boundary_info
# gradients of shape functions
if g.dim == 2:
q, w = g.reference_element.quadrature(order=2)
SF_GRAD = np.array(([q[..., 1] - 1., q[..., 0] - 1.],
[1. - q[..., 1], -q[..., 0]],
[q[..., 1], q[..., 0]],
[-q[..., 1], 1. - q[..., 0]]))
else:
raise NotImplementedError
q, w = g.reference_element.quadrature(order=2)
SF_GRAD = LagrangeShapeFunctionsGrads[g.reference_element][1]
SF_GRAD = SF_GRAD(q)
self.logger.info('Calulate gradients of shape functions transformed by reference map ...')
self.logger.info('Calculate gradients of shape functions transformed by reference map ...')
SF_GRADS = np.einsum('eij,pjc->epic', g.jacobian_inverse_transposed(0), SF_GRAD)
self.logger.info('Calculate all local scalar products beween gradients ...')
self.logger.info('Calculate all local scalar products between gradients ...')
if self.diffusion_function is not None and self.diffusion_function.shape_range == ():
D = self.diffusion_function(self.grid.centers(0), mu=mu)
SF_INTS = np.einsum('epic,eqic,c,e,e->epq', SF_GRADS, SF_GRADS, w, g.integration_elements(0), D).ravel()
......@@ -652,7 +633,7 @@ class AdvectionOperatorP1(NumpyMatrixBasedOperator):
(Lu)(x) = c ∇ ⋅ [ v(x) u(x) ]
The function `v` is vector-valued.
The function `v` has to be vector-valued.
Parameters
----------
......@@ -695,28 +676,18 @@ class AdvectionOperatorP1(NumpyMatrixBasedOperator):
g = self.grid
bi = self.boundary_info
# gradients of shape functions
if g.dim == 2:
SF_GRAD = np.array(([-1., -1.],
[1., 0.],
[0., 1.]))
SF = [lambda X: 1 - X[..., 0] - X[..., 1],
lambda X: X[..., 0],
lambda X: X[..., 1]]
# SF_GRAD(function, component)
else:
raise NotImplementedError
q, w = g.reference_element.quadrature(order=1)
SF = LagrangeShapeFunctions[g.reference_element][1]
self.logger.info('Calulate gradients of shape functions transformed by reference map ...')
self.logger.info('Calculate gradients of shape functions transformed by reference map ...')
SF_GRAD = LagrangeShapeFunctionsGrads[g.reference_element][1]
SF_GRADS = np.einsum('eij,pj->epi', g.jacobian_inverse_transposed(0), SF_GRAD)
# SF_GRADS(element, function, component)
SFQ = np.array(tuple(f(q) for f in SF))
# SFQ(function, quadraturepoint)
self.logger.info('Calculate all local scalar products beween gradients ...')
self.logger.info('Calculate all local scalar products between gradients ...')
D = self.advection_function(self.grid.centers(0), mu=mu)
SF_INTS = - np.einsum('pc,eqi,c,e,ei->eqp', SFQ, SF_GRADS, w, g.integration_elements(0), D).ravel()
del D
......@@ -804,29 +775,17 @@ class AdvectionOperatorQ1(NumpyMatrixBasedOperator):
g = self.grid
bi = self.boundary_info
# gradients of shape functions
if g.dim == 2:
q, w = g.reference_element.quadrature(order=2)
SF_GRAD = np.array(([q[..., 1] - 1., q[..., 0] - 1.],
[1. - q[..., 1], -q[..., 0]],
[q[..., 1], q[..., 0]],
[-q[..., 1], 1. - q[..., 0]]))
# SF_GRAD(function, component, quadraturepoint)
SF = [lambda X: (1 - X[..., 0]) * (1 - X[..., 1]),
lambda X: (1 - X[..., 1]) * (X[..., 0]),
lambda X: (X[..., 0]) * (X[..., 1]),
lambda X: (1 - X[..., 0]) * (X[..., 1])]
else:
raise NotImplementedError
self.logger.info('Calulate gradients of shape functions transformed by reference map ...')
self.logger.info('Calculate gradients of shape functions transformed by reference map ...')
q, w = g.reference_element.quadrature(order=2)
SF_GRAD = LagrangeShapeFunctionsGrads[g.reference_element][1](q)
SF_GRADS = np.einsum('eij,pjc->epic', g.jacobian_inverse_transposed(0), SF_GRAD)
# SF_GRADS(element,function,component,quadraturepoint)
SF = LagrangeShapeFunctions[g.reference_element][1]
SFQ = np.array(tuple(f(q) for f in SF))
# SFQ(function, quadraturepoint)
self.logger.info('Calculate all local scalar products beween gradients ...')
self.logger.info('Calculate all local scalar products between gradients ...')
D = self.advection_function(self.grid.centers(0), mu=mu)
SF_INTS = - np.einsum('pc,eqic,c,e,ei->eqp', SFQ, SF_GRADS, w, g.integration_elements(0), D).ravel()
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment