Skip to content
Snippets Groups Projects
Commit 02b2a227 authored by Laura Christine Kühle's avatar Laura Christine Kühle
Browse files

Moved 'do_initial_projection()' to projection_utils.py.

parent 2c7ae512
Branches
No related tags found
No related merge requests found
...@@ -7,8 +7,8 @@ import os ...@@ -7,8 +7,8 @@ import os
import time import time
import numpy as np import numpy as np
from .DG_Approximation import do_initial_projection
from .Mesh import Mesh from .Mesh import Mesh
from .projection_utils import do_initial_projection
from .Quadrature import Gauss from .Quadrature import Gauss
from .Basis_Function import OrthonormalLegendre from .Basis_Function import OrthonormalLegendre
......
...@@ -4,20 +4,17 @@ ...@@ -4,20 +4,17 @@
""" """
import json import json
import numpy as np
from sympy import Symbol
import math import math
from .Basis_Function import Basis from .Basis_Function import Basis
from .encoding_utils import encode_ndarray from .encoding_utils import encode_ndarray
from .Initial_Condition import InitialCondition from .Initial_Condition import InitialCondition
from .Mesh import Mesh from .Mesh import Mesh
from .projection_utils import do_initial_projection
from .Quadrature import Quadrature from .Quadrature import Quadrature
from .Troubled_Cell_Detector import TroubledCellDetector from .Troubled_Cell_Detector import TroubledCellDetector
from .Update_Scheme import UpdateScheme from .Update_Scheme import UpdateScheme
x = Symbol('x')
class DGScheme: class DGScheme:
"""Class for Discontinuous Galerkin Method. """Class for Discontinuous Galerkin Method.
...@@ -152,59 +149,3 @@ class DGScheme: ...@@ -152,59 +149,3 @@ class DGScheme:
with open(data_file + '.json', 'w') \ with open(data_file + '.json', 'w') \
as json_file: as json_file:
json_file.write(json.dumps(approx_stats)) json_file.write(json.dumps(approx_stats))
def do_initial_projection(init_cond, mesh, basis, quadrature,
x_shift=0):
"""Calculates initial projection.
Calculates a projection at time step 0 and adds ghost cells on both
sides of the array.
Parameters
----------
init_cond : InitialCondition object
Initial condition used for calculation.
mesh : Mesh object
Mesh for calculation.
basis: Basis object
Basis used for calculation.
quadrature: Quadrature object
Quadrature used for evaluation.
x_shift: float, optional
Shift in x-direction for each evaluation point. Default: 0.
Returns
-------
ndarray
Matrix containing projection of size (N+2, p+1) with N being the
number of mesh cells and p being the polynomial degree of the basis.
"""
# Initialize output matrix
output_matrix = np.zeros((basis.polynomial_degree+1,
mesh.num_cells+2*mesh.num_ghost_cells))
# Calculate basis based on quadrature
basis_matrix = np.array([np.vectorize(basis.basis[degree].subs)(
x, quadrature.nodes) for degree in range(basis.polynomial_degree+1)])
# Calculate points based on initial condition
points = quadrature.nodes[:, np.newaxis] @ \
np.repeat(mesh.cell_len/2, mesh.num_cells)[:, np.newaxis].T + \
np.tile(mesh.non_ghost_cells - x_shift, (quadrature.num_nodes, 1))
init_matrix = np.vectorize(init_cond.calculate, otypes=[np.float])(
x=points, mesh=mesh)
# Set output matrix for regular cells
output_matrix[:, mesh.num_ghost_cells:
output_matrix.shape[-1]-mesh.num_ghost_cells] = \
(basis_matrix * quadrature.weights) @ init_matrix
# Set ghost cells
output_matrix[:, :mesh.num_ghost_cells] = \
output_matrix[:, -2*mesh.num_ghost_cells:-mesh.num_ghost_cells]
output_matrix[:, output_matrix.shape[-1]-mesh.num_ghost_cells:] = \
output_matrix[:, mesh.num_ghost_cells:2*mesh.num_ghost_cells]
return output_matrix
...@@ -91,3 +91,59 @@ def calculate_exact_solution( ...@@ -91,3 +91,59 @@ def calculate_exact_solution(
exact = np.reshape(exact, (1, exact.size)) exact = np.reshape(exact, (1, exact.size))
return grid, exact return grid, exact
def do_initial_projection(init_cond, mesh, basis, quadrature,
x_shift=0):
"""Calculates initial projection.
Calculates a projection at time step 0 and adds ghost cells on both
sides of the array.
Parameters
----------
init_cond : InitialCondition object
Initial condition used for calculation.
mesh : Mesh object
Mesh for calculation.
basis: Basis object
Basis used for calculation.
quadrature: Quadrature object
Quadrature used for evaluation.
x_shift: float, optional
Shift in x-direction for each evaluation point. Default: 0.
Returns
-------
ndarray
Matrix containing projection of size (N+2, p+1) with N being the
number of mesh cells and p being the polynomial degree of the basis.
"""
# Initialize output matrix
output_matrix = np.zeros((basis.polynomial_degree+1,
mesh.num_cells+2*mesh.num_ghost_cells))
# Calculate basis based on quadrature
basis_matrix = np.array([np.vectorize(basis.basis[degree].subs)(
x, quadrature.nodes) for degree in range(basis.polynomial_degree+1)])
# Calculate points based on initial condition
points = quadrature.nodes[:, np.newaxis] @ \
np.repeat(mesh.cell_len/2, mesh.num_cells)[:, np.newaxis].T + \
np.tile(mesh.non_ghost_cells - x_shift, (quadrature.num_nodes, 1))
init_matrix = np.vectorize(init_cond.calculate, otypes=[np.float])(
x=points, mesh=mesh)
# Set output matrix for regular cells
output_matrix[:, mesh.num_ghost_cells:
output_matrix.shape[-1]-mesh.num_ghost_cells] = \
(basis_matrix * quadrature.weights) @ init_matrix
# Set ghost cells
output_matrix[:, :mesh.num_ghost_cells] = \
output_matrix[:, -2*mesh.num_ghost_cells:-mesh.num_ghost_cells]
output_matrix[:, output_matrix.shape[-1]-mesh.num_ghost_cells:] = \
output_matrix[:, mesh.num_ghost_cells:2*mesh.num_ghost_cells]
return output_matrix
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment