Commit f5fe1382 authored by Élise Le Mélédo's avatar Élise Le Mélédo
Browse files

Added new quadrature rules, updated docs, added safe mesh runtime export

parent 38356879
......@@ -2,6 +2,10 @@
# (its number of edges basically), the order and the barycentric
# coordinates of the considered point with respect to the given triangle
# ==============================================================================
# Pythonic interfacing constraints (safety barrier)
# ==============================================================================
__all__ = ["BarycentricBasis","BarycentricGradients"]
# ==============================================================================
# Preliminary imports
......@@ -10,7 +14,60 @@ import numpy as np
import itertools
# ==============================================================================
# Basis functions for triangular shape
# Interfacing the basis functions depending on the element's type
# ==============================================================================
#### Interfacing routine for the basis function computations ####
def BarycentricBasis(Nele, Type, Order, points):
"""
Simple mapping to the implementation giving out the barycentic basis functions according to the elment's type (simplicial, quad, ...)
Args:
Nele (integer): the number of edges of the considered element to compute the basis functions for
Type (string): type of wished basis functions
Order (integer): order of the wished basis functions
points (2D numpy array): the points at which to evaluate the basis functons (NbPoints, 3)
Returns:
EvaluatedBase (2D numpy array): the value of all the basis functions at each given point (NbBasisFunc x NbPoints)
through the call of the relevant function.
"""
# ---------------- Mapping to the right function -------------------------------------------------
# Case of a simplicial shape
if Nele==3: return(SimplicialBasis(Type, Order, points))
else: raise(ValueError("Error: Barycentric basis functions are not yet implemented for elements having "+str(Nele)+"edges. Abortion."))
#### Interfacing routine for the gradient computations ####
def BarycentricGradients(Nele, Type, Order, points, n, vol):
"""
Simple mapping to the implementation giving out the barycentic basis functions gradients according to the elment's type (simplicial, quad, ...)
Args:
Nele (integer): the number of edges of the considered element to compute the basis functions for
Type (string): type of wished basis functions
Order (integer): order of the wished basis functions
points (2D numpy array): the points at which to evaluate the basis functons (NbPoints, 3)
n (2D numpy array): the normal vector to each face, in the order corresponding to the vertices furnishing the barycentric coordinates referential
vol (float): the volume of the element furnishing the barycentric coordinates referential
Returns:
EvaluatedGradient (3D numpy array): the value of all the x and y gradient's coordinates at each given point, ordered accordingly to the
vertices from which the barycentic coordinates have been generated (NbBasisFunc x NbPoints x 2),
through the call of the relevant function.
"""
# ---------------- Mapping to the right function -------------------------------------------------
# Case of a simplicial shape
if Nele==3: return(SimplicialGradients(Type, Order, points, n, vol))
else: raise(ValueError("Error: Barycentric basis functions are not yet implemented for elements having "+str(Nele)+"edges. Abortion."))
################################################################################
################################################################################
# ==============================================================================
# Technical routines for the basis functions computations
# ==============================================================================
#### Simplicial basis function from barycentric coordinates ####
......@@ -53,6 +110,11 @@ def SimplicialBasis(Type, Order, points):
# ------------------ Return the considered basis ----------------------------------
return(EvaluatedBase)
# ==============================================================================
# Technical routines for the basis functions gradient computations
# ==============================================================================
#### Simplicial gradient of basis function from barycentric coordinates ####
def SimplicialGradients(Type, Order, points, n, vol):
"""
......@@ -101,7 +163,7 @@ def SimplicialGradients(Type, Order, points, n, vol):
# Scale the values to the geometry
Grad = np.dot(np.swapaxes(np.array([EvaluatedGradient[:,:,0]-EvaluatedGradient[:,:,2],
EvaluatedGradient[:,:,1]-EvaluatedGradient[:,:,2]]),0,2),\
n[1:,:])/(2.*vol)
n[1:,:])/(2.*vol)
# Safety check
else: raise NotImplementedError("The type of wished basis function is not implemented yet. Abortion.")
......
# File that procures the necessary tools to handle barycentric coordinates
# File that procures the necessary tools to handle barycentric coordinates
# ==============================================================================
......@@ -13,30 +13,29 @@ import numpy as np
#### Get the barycentric coordinates in their classical definiton ####
def GetBarycentricCoordinates(points, vertices):
""" Barycentric coordinates from the given vertices (forming a triangles)
Args:
Points (2D numpy array): Cartesian coordinates of the points to convert (NbPoints x 2)
vertices (2D numpy array): Cartesian coordinates of the points taken as reference, ordered. (3 x 2)
Returns:
bcoors (2D numpy array): Barycentric coordinates of points with respect to the vertices (NbPoints x 3)
----------------------------------------------------------------------"""
"""
# Performing the computations of the barycentric coordinates
hyp1 = vertices[1,:]-vertices[0,:]
hyp2 = vertices[2,:]-vertices[0,:]
dist = points-vertices[0,:]
dhyp11 = np.dot(hyp1, hyp1)
dhyp12 = np.dot(hyp1, hyp2)
dhyp22 = np.dot(hyp2, hyp2)
dhyp31 = np.dot(dist, hyp1)
dhyp32 = np.dot(dist, hyp2)
quotu = dhyp11*dhyp22-dhyp12**2
b1 = (dhyp22*dhyp31-dhyp12*dhyp32)/quotu
b2 = (dhyp11*dhyp32-dhyp12*dhyp31)/quotu
bcoords = np.array([1-b1-b2, b1, b2]).T
......@@ -47,17 +46,15 @@ def GetBarycentricCoordinates(points, vertices):
#### Retrieving the cartesian coordinates ####
def GetCartesianCoordinates(bcoords, vertices):
""" Barycentric coordinates from the given vertices (forming a triangles)
Args:
bcoors (2D numpy array): Barycentric coordinates of points with respect to the vertices (NbPoints x 3)
vertices (2D numpy array): Cartesian coordinates of the points taken as reference, ordered. (3 x 2)
Returns:
Points (2D numpy array): Cartesian coordinates of the points to convert (NbPoints x 2)
----------------------------------------------------------------------"""
"""
# Computing back the cartesian points and returning them
points = np.dot(bcoords, vertices)
return(points)
#Module furnishing the lobatto quadrature points and weights, copied
#from an external library belonging to Greg von Winckel - 04/17/2004
#from an external library belonging to Greg von Winckel - 04/17/2004
# ==============================================================================
......@@ -15,7 +15,7 @@ def GaussLobatto(n,eps):
Computes the Legendre-Gauss-Lobatto nodes, weights and the LGL Vandermonde
matrix. The LGL nodes are the zeros of (1-x^2)*P'_N(x). Useful for numerical
integration and spectral methods.
Initial content from Greg von Winckel - 04/17/2004, Python translation of lglnodes.m,
Initial content from Greg von Winckel - 04/17/2004, Python translation of lglnodes.m,
translated and modified into Python by Jacob Schroder - 9/15/2018
......@@ -24,27 +24,25 @@ def GaussLobatto(n,eps):
epse (float): safety net for division
Returns:
(nodes, weights) (tuple): representing the quadrature nodes and weights.
(nodes, weights) (tuple): representing the quadrature nodes and weights.
.. rubric:: Example
.. code::
.. code::
>>> from lglnodes import *
>>> (nodes, weights) = lglnodes(3)
>>> print(str(nodes) + " " + str(weights))
[-1. -0.4472136 0.4472136 1. ] [0.16666667 0.83333333 0.83333333 0.16666667]
.. note::
- (n+1) nodes and weights are returned
- Reference on LGL nodes and weights: C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Tang, "Spectral Methods in Fluid Dynamics," Section 2.3. Springer-Verlag 1987
"""
# Initialisation
# -------------------- Initialisation --------------------------------------
w = np.zeros((n+1,))
x = np.zeros((n+1,))
xold = np.zeros((n+1,))
......@@ -52,6 +50,7 @@ def GaussLobatto(n,eps):
# The Legendre Vandermonde Matrix
P = np.zeros((n+1,n+1))
# -------------------- Computations ----------------------------------------
# Use the Chebyshev-Gauss-Lobatto nodes as the first guess
for i in range(n+1): x[i] = -np.cos(np.pi*i / n)
......@@ -69,6 +68,9 @@ def GaussLobatto(n,eps):
if (np.max(np.abs(x - xold).flatten()) < eps): break
# Compute the weights
w = 2.0/((n*(n+1))*(P[:,n]**2))
# -------------------- End-tasks -------------------------------------------
# Return the values
return(x, w)
# File that procures simple pythonic hacks
# ==============================================================================
# Pythonic hack routines
# ==============================================================================
#### Partial function applying the last arguments ####
def partial(func, *args):
"""
Partial function applying the last arguments
Partial function applying the last arguments
Args:
func (function callback): the function of interest
*args (arguments list): the arguments to parse to the function from the end.
Returns:
Lmb (function callback): the lambda function defining the partial callback from the last arguments
"""
return(lambda *a: func(*(a + args)))
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
################################################################################
####################################################################################
## ##
## Problem describing a simple constant rotation of a bump in a given fluid. ##
## ##
####################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
......@@ -28,7 +31,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = [] #(2,0.6,1.5,0,0.5),(4,0.5, 1.0, 0.8, 0.8, 0.05, 0.15)] #(4,0.5, 1.0, 0.8, 0.5, 0.05, 0.15) [[[0.5, 0],[0.5, 0.5],[0, 0.5]],[[0.75,0.75],[0.75,1],[1,1]], 2, (2,3), (1,3,5,6)] (0,9,0,0.5,-0.5,0.5,0.5)
self.FluidInitialSubdomain = []
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(2,)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
################################################################################
####################################################################################
## ##
## Problem describing a simple constant translation of a bump in a given fluid. ##
## ##
####################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
......@@ -28,7 +31,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = [] #(2,0.6,1.5,0,0.5),(4,0.5, 1.0, 0.8, 0.8, 0.05, 0.15)] #(4,0.5, 1.0, 0.8, 0.5, 0.05, 0.15) [[[0.5, 0],[0.5, 0.5],[0, 0.5]],[[0.75,0.75],[0.75,1],[1,1]], 2, (2,3), (1,3,5,6)] (0,9,0,0.5,-0.5,0.5,0.5)
self.FluidInitialSubdomain = []
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(1,)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
## ##
## Problem describing a simple bump advection with two coexisting fluids ##
## propagating at different speed to see the behaviour of the level set when ##
## the bump of fluid one enters the fluid 2. No particular physical meaning. ##
## ##
################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
# ******************************************************************************
# Meshing
#*******************************************************************************
self.ProblemID = "Advection_Translation_2Fluids" # Should match the module name
self.ProblemID = "Advection_Translation_2Fluids"
self.MeshName = "BigSimpleRectangle"
# ******************************************************************************
......@@ -28,7 +33,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = [5] #(2,0.6,1.5,0,0.5),(4,0.5, 1.0, 0.8, 0.8, 0.05, 0.15)] #(4,0.5, 1.0, 0.8, 0.5, 0.05, 0.15) [[[0.5, 0],[0.5, 0.5],[0, 0.5]],[[0.75,0.75],[0.75,1],[1,1]], 2, (2,3), (1,3,5,6)] (0,9,0,0.5,-0.5,0.5,0.5)
self.FluidInitialSubdomain = [5]
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(0,0.1), (1,)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
## ##
## Problem describing a simple contant advection with two coexisting fluids ##
## having the same state (height). Only for simple debug test, no physical ##
## meaning. ##
## ##
################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
# ******************************************************************************
# Meshing
#*******************************************************************************
self.ProblemID = "Advection_Translation_Cste_2Fluids" # Should match the module name
self.ProblemID = "Advection_Translation_Cste_2Fluids"
self.MeshName = "BigSimpleSquare"
# ******************************************************************************
......@@ -28,7 +33,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = [5] #(2,0.6,1.5,0,0.5),(4,0.5, 1.0, 0.8, 0.8, 0.05, 0.15)] #(4,0.5, 1.0, 0.8, 0.5, 0.05, 0.15) [[[0.5, 0],[0.5, 0.5],[0, 0.5]],[[0.75,0.75],[0.75,1],[1,1]], 2, (2,3), (1,3,5,6)] (0,9,0,0.5,-0.5,0.5,0.5)
self.FluidInitialSubdomain = [5]
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(0,1), (0,1)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
## ##
## Problem describing a simple translation of a bump evolving in a fluid 1 ##
## towards a fluid2. The associated mesh contains a hole to perturb the ##
## behaviour. ##
## ##
################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
# ******************************************************************************
# Meshing
#*******************************************************************************
self.ProblemID = "Advection_Translation_Obstacle_2Fluids" # Should match the module name
self.ProblemID = "Advection_Translation_Obstacle_2Fluids"
self.MeshName = "LShapeWithHoleAndRoundCorner"
# ******************************************************************************
......@@ -29,7 +34,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = [(4, 0.18, 0.5, 0.14, 0.14, 0.14, 1, 1)]
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(0,0), (0,0.5)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
## ##
## Problem describing a the Karni problem with a bubble of Helium ##
## ##
################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
......@@ -29,7 +32,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = [(4,0,0,0.5,0.5,0.5,4,4)]
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(7, 1.22), (6,)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
## ##
## Problem describing a the Karni problem with a bubble of R22 ##
## ##
################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
......@@ -29,7 +32,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = [(4,0,0,0.5,0.5,0.5,4,4)]
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(7, 1.22), (6,)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
## ##
## Problem describing a the Karni problem with a bubble of Helium whose ##
## shape is artificially set to a square. ##
## ##
################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
......@@ -29,7 +33,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = [[[-0.25, -0.250],[0.25, -0.25],[0.25, 0.25],[-0.25,0.25]]]
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(7, 1.22), (6,)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
## ##
## Problem describing a SOD problem with a unique fluid in a big square ##
## domain. ##
## ##
################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
# ******************************************************************************
# Meshing
#*******************************************************************************
self.ProblemID = "Euler_SOD" # Should match the module name
self.ProblemID = "Euler_SOD"
self.MeshName = "SimpleSquare"
# ******************************************************************************
......@@ -29,7 +33,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = []
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(3,)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
## ##
## Problem describing a SOD problem with a one fluid type per initial ##
## submdomain. ##
## ##
################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
......@@ -29,7 +33,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = [(4,0,0,0.5,0.5,0.5,4,4)]
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(5,), (4,)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
## ##
## Problem describing a SOD problem with a unique fluid in a small ##
## rectangular domain. ##
## ##
################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
......@@ -29,7 +33,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = [(4,0,0,0.5,0.5,0.5,4,4)]
# Index of the initial conditions to apply on each subdomain
self.InitialConditions = [(5,), (4,)]
......
################################################################################
##
## Template file to design the properties of the problem of interest ##
##
## ##
## Problem describing a Shock tube problem with a Mach number of 1.22 ##
## (with Euler governing equations) ##
## ##
################################################################################
# ------------- Preliminary imports --------------------
import numpy as np
# ------------- Problem definition ---------------------
class ProblemData():
def __init__(self):
......@@ -29,7 +33,7 @@ class ProblemData():
# The fluid at index 0 is the complementary to the specified ones.
# Only need to precise the location of other fluids
self.FluidInitialSubdomain = []
# Index of the initial conditions to apply on each subdomain