Skip to content
Snippets Groups Projects
Commit 9693c9d3 authored by Edilbert Christhuraj's avatar Edilbert Christhuraj
Browse files

Initial commit message

parents
Branches
Tags WS2324
No related merge requests found
Pipeline #412442 skipped
Stage: build
Hide whitespace changes
Inline Side-by-side
Showing
with 1346 additions and 0 deletions
.dockerignore 0 → 100644
+ 1
− 0
View file @ 9693c9d3
**/.git
.gitignore 0 → 100644
+ 16
− 0
View file @ 9693c9d3
#ignore FILE TYPES
#*.xml #not relevant
#*.msh #not relevant
*.h5
*.xdmf
#ignore IPython
*.ipynb_checkpoints
*__pycache__*
#ignore FOLDERS
etc/
results_mathematica/
input/
thermal_edge_flow/
###############################################
.gitlab-ci.yml 0 → 100644
+ 25
− 0
View file @ 9693c9d3
image: docker:stable
variables:
DOCKER_TLS_CERTDIR: ""
DOCKER_HOST: tcp://docker:2375
DOCKER_DRIVER: overlay2
services:
- docker:dind # docker in docker
stages:
- build
# **************************************************************************** #
# build
# **************************************************************************** #
build:environment:
stage: build
before_script:
- docker login -u $CI_REGISTRY_USER -p $CI_REGISTRY_PASSWORD $CI_REGISTRY
script:
- docker pull $CI_REGISTRY_IMAGE:latest || true
- docker build --cache-from $CI_REGISTRY_IMAGE:latest --tag $CI_REGISTRY_IMAGE:$CI_COMMIT_SHA --tag $CI_REGISTRY_IMAGE:latest .
# - docker push $CI_REGISTRY_IMAGE:$CI_COMMIT_SHA # no need
- docker push $CI_REGISTRY_IMAGE:latest
when: manual
CHANGELOG.rst 0 → 100644
+ 2
− 0
View file @ 9693c9d3
Describe the state
\ No newline at end of file
CODEOWNERS 0 → 100644
+ 3
− 0
View file @ 9693c9d3
# default match
* @19ec94
\ No newline at end of file
CONTRIBUTING.rst 0 → 100644
+ 2
− 0
View file @ 9693c9d3
Mention the contributors
\ No newline at end of file
Dockerfile 0 → 100644
+ 37
− 0
View file @ 9693c9d3
# Use FEniCS base image
FROM quay.io/fenicsproject/stable:2019.1.0.r3
# Descriptions
LABEL maintainer="Edilbert Christhuraj <bert.edil@gmail.com>"
LABEL description="FEnicS-Rx"
# Specify software versions
ENV GMSH_VERSION 4.4.0
# Download Install Gmsh SDK with dependecies from Github's dolfinx Dockerfile
RUN export DEBIAN_FRONTEND=noninteractive && \
apt-get -qq update && \
# apt-get -yq --with-new-pkgs -o Dpkg::Options::="--force-confold" upgrade && \ # skip 110 package updates/upgrades
apt-get -y install \
libglu1 \
libxcursor-dev \
libxinerama1 && \
apt-get clean && \
rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/*
RUN cd /usr/local && \
wget -nc http://gmsh.info/bin/Linux/gmsh-${GMSH_VERSION}-Linux64-sdk.tgz && \
tar -xf gmsh-${GMSH_VERSION}-Linux64-sdk.tgz
ENV PATH=/usr/local/gmsh-${GMSH_VERSION}-Linux64-sdk/bin:$PATH
# Install any needed packages specified in requirements.txt
# RUN pip install --trusted-host pypi.python.org -r requirements.txt
COPY requirements.txt /tmp/
RUN pip3 install --requirement /tmp/requirements.txt
# Install the fenicsR13 package (puts it into the PATH)
#COPY . /tmp/
#RUN pip install --editable /tmp/.
# Replace default FEniCS Docker WELCOME screen with custom WELCOME screen
COPY WELCOME .
RUN echo "Built: $(date)" >> WELCOME
LICENSE 0 → 100644
+ 2
− 0
View file @ 9693c9d3
Add license
\ No newline at end of file
README.rst 0 → 100644
+ 11
− 0
View file @ 9693c9d3
FFME: A generic framework to solve arbitrary set of moment equations arise in the context of moment methods in kinetic gas theory
=================================================================================================================================
Main features
-------------------------------------------------------------------------------
-blabla
-blabla blabla
Installation
-------------
WELCOME 0 → 100644
+ 10
− 0
View file @ 9693c9d3
--------------------------------------------------------------------------------
-> Contact: Edilbert Christhuraj <Christhuraj@mathcces.rwth-aachen.de>
-> Contact: Manuel Torrilhon <mt@mathcces.rwth-aachen.de>
-> Repository: <https://git.rwth-aachen.de/bert.edil/fenics_rx>
Welcome to FEniCS-Rx
The Example files can be found in /shared/examples/
--------------------------------------------------------------------------------
docker-compose.yml 0 → 100644
+ 36
− 0
View file @ 9693c9d3
version: '3.2'
services:
fenics_rx_from_git:
image: registry.git.rwth-aachen.de/bert.edil/fenics_rx:latest
working_dir: /home/fenics/shared
volumes:
- type: bind
source: .
target: /home/fenics/shared
- type: volume
source: instant-cache
target: /home/fenics/.cache
volume:
nocopy: true
stdin_open: true
tty: true
fenics_rx_local:
build:
context: .
dockerfile: Dockerfile
working_dir: /home/fenics/shared
volumes:
- type: bind
source: .
target: /home/fenics/shared
- type: volume
source: instant-cache
target: /home/fenics/.cache
volume:
nocopy: true
stdin_open: true
tty: true
volumes:
instant-cache:
examples/ffme/ffme.ipynb 0 → 100644
+ 425
− 0
View file @ 9693c9d3
%% Cell type:code id: tags:
``` python
import dolfin as df
import numpy as np
import os
import time as time_module
import matplotlib.pyplot as plt
import csv
import yaml #version higher than 5.0
df.parameters['ghost_mode']= 'shared_facet' #for parallel program
#My library
from modules.SystemMatrix import SystemMatrix # <--System Matrix class
from modules.Mesh import H5Mesh # <-- Mesh class
from modules.FunctionSpace import FunctionSpace #<--Function Space Class
from modules.FormulateVariationalProblem import FormulateVariationalProblem
from modules.Solver import Solver
#Inputs
with open(r'input_r13_ring_nosource_rot_inflow.yml') as input_file:
my_input_cls = yaml.load(input_file, Loader=yaml.FullLoader)
projection = False
manual = True
automatic = False
my_current_mesh = my_input_cls['mesh_list'][0]
#step-1 : Convert current_mesh_name into a char_list
#step-2: Select the 4th-to-last element in the char_list
#Step-3: Convert the 4th-to-last element to integer
#Now I can run any single mesh I want in whichever order
mesh_number=0 #For saving files according to mesh number
mesh_number = int(list(my_current_mesh)[-4])
print(mesh_number)
```
%% Output
1
%% Cell type:code id: tags:
``` python
#MESH
my_mesh_cls = H5Mesh(my_current_mesh) #<-- Class instance is created
#print(my_current_mesh)
print("Maximum value of the mesh: ", my_mesh_cls.mesh.hmax())
df.plot(my_mesh_cls.mesh)
```
%% Output
Maximum value of the mesh: 0.6340332990696501
[<matplotlib.lines.Line2D at 0x7f4d629ea978>,
<matplotlib.lines.Line2D at 0x7f4d629eab00>]
%% Cell type:code id: tags:
``` python
#FUNCTION SPACE
problem_type = my_input_cls['problem_type']
number_of_moments = my_input_cls['number_of_moments']
my_function_space_cls = FunctionSpace(my_input_cls, my_mesh_cls) #<-- class instance
my_function_space_cls.set_function_space()
```
%% Cell type:code id: tags:
``` python
#SYSTEM MATRIX
my_system_matrix_cls = SystemMatrix(my_input_cls, my_mesh_cls) # <-- class instance
my_system_matrix_cls.convert_to_ufl_form()
```
%% Cell type:code id: tags:
``` python
#my_system_matrix_cls.SP_Coeff
```
%% Cell type:code id: tags:
``` python
#my_system_matrix_cls.SA_x
#np.savetxt("SAx_old.txt", my_system_matrix_cls.SA_x, fmt="%d")
#np.savetxt("SAy_old.txt", my_system_matrix_cls.SA_y, fmt="%d")
#np.savetxt("SP_Coeff_old.txt", my_system_matrix_cls.SP_Coeff, fmt="%d")
#np.savetxt("BC1_old.txt", my_system_matrix_cls.BC1, fmt="%s")
#np.savetxt("BC2_old.txt", my_system_matrix_cls.BC2, fmt="%s")
#np.savetxt("BC1_rhs_old.txt", my_system_matrix_cls.BC1_rhs, fmt="%s")
#np.savetxt("BC2_rhs_old.txt", my_system_matrix_cls.BC2_rhs, fmt="%s")
```
%% Cell type:code id: tags:
``` python
#VARIATIONAL PROBLEM
my_var_prob_cls = FormulateVariationalProblem(
my_input_cls,
my_mesh_cls,
my_system_matrix_cls,
my_function_space_cls
) #<-- class instance
my_var_prob_cls.create_lhs()
my_var_prob_cls.create_rhs()
```
%% Cell type:code id: tags:
``` python
#SOLVER
my_solver_cls = Solver(my_input_cls, my_function_space_cls, my_var_prob_cls) #<-- class instance
if problem_type == 'nonlinear':
my_solver_cls.custom_newton_solver()
u = my_solver_cls.u
else:
my_solver_cls.inbuilt_linear_solver()
u = my_solver_cls.u_Function
```
%% Cell type:code id: tags:
``` python
#===============
#Post-processing
#===============
sol = u.split()
```
%% Cell type:code id: tags:
``` python
c =df.plot(sol[3])
plt.colorbar(c)
plt.title('Temparature')
plt.savefig('edge.png', bbox_inches='tight')
```
%% Output
%% Cell type:code id: tags:
``` python
def write_func(field_name, variable_name, mesh_num):
xdmffile_u = df.XDMFFile(df.MPI.comm_world,
'results_mathematica/{0}_{1}.xdmf'
.format(variable_name,mesh_num))
xdmffile_u.write(field_name)
xdmffile_u.close()
for i in range(my_input_cls["number_of_moments"]):
write_func(sol[i], i, mesh_number)
```
%% Cell type:code id: tags:
``` python
#file = File("sol.pvd")
#file << u.split()[0]
#file << u.split()[1]
#file << u.split()[2]
'''
file = df.File("sol0.pvd")
file.write(sol[0])
file1 = df.File("sol1.pvd")
file1.write(sol[1])
file2 = df.File("sol2.pvd")
file2.write(sol[2])
'''
```
%% Output
'\nfile = df.File("sol0.pvd")\nfile.write(sol[0])\nfile1 = df.File("sol1.pvd")\nfile1.write(sol[1])\nfile2 = df.File("sol2.pvd")\nfile2.write(sol[2])\n'
%% Cell type:code id: tags:
``` python
xmesh = my_function_space_cls.V.tabulate_dof_coordinates()[:,0]
ymesh = my_function_space_cls.V.tabulate_dof_coordinates()[:,1]
ux =sol[1].vector().get_local()
uy =sol[2].vector().get_local()
```
%% Cell type:code id: tags:
``` python
plt.quiver(xmesh,ymesh,ux,uy, color='Teal', headlength=1)
```
%% Output
<matplotlib.quiver.Quiver at 0x7f4d5fe6bd68>
%% Cell type:code id: tags:
``` python
#X_, Y_ = np.meshgrid(xmesh,ymesh)
'''
#Z = []
for i in range(len(xmesh)):
new =[]
for i in range(len(xmesh)):
new.append([sol[1].vector().get_local()[i], sol[2].vector().get_local()[j]])
Z.append(new)
Z = np.array(Z)
'''
```
%% Output
'\n#Z = []\nfor i in range(len(xmesh)):\n new =[]\n for i in range(len(xmesh)):\n new.append([sol[1].vector().get_local()[i], sol[2].vector().get_local()[j]])\n Z.append(new)\nZ = np.array(Z)\n'
%% Cell type:code id: tags:
``` python
plt.tripcolor(xmesh,ymesh,ux)
#plt.tricontour(xmesh,ymesh,uy)
```
%% Output
<matplotlib.collections.PolyCollection at 0x7f4d5fc719e8>
%% Cell type:code id: tags:
``` python
if automatic:
u_old = u
u_old.split()[3].compute_vertex_values()-u.split()[3].compute_vertex_values()
```
%% Output
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-17-791f8ef241af> in <module>()
1 if automatic:
2 u_old = u
----> 3 u_old.split()[3].compute_vertex_values()-u.split()[3].compute_vertex_values()
NameError: name 'u_old' is not defined
examples/ffme/ffme.py 0 → 100644
+ 213
− 0
View file @ 9693c9d3
import sys
import dolfin as df
import numpy as np
import os
import time as time_module
import matplotlib.pyplot as plt
import csv
import yaml #version higher than 5.0
df.parameters['ghost_mode']= 'shared_facet' #for parallel program
#My library
from modules.SystemMatrix import SystemMatrix # <--System Matrix class
from modules.Mesh import H5Mesh # <-- Mesh class
from modules.FunctionSpace import FunctionSpace #<--Function Space Class
from modules.FormulateVariationalProblem import FormulateVariationalProblem
from modules.Solver import Solver
#Inputs
if len(sys.argv) >1 :
user_given_input_file= sys.argv[1]
else:
print("Provide an input file")
quit()
with open(user_given_input_file,'r') as input_file:
my_input_cls = yaml.load(input_file, Loader=yaml.FullLoader)
#projection = False
#manual = False
#automatic = True
for current_mesh in range(len(my_input_cls["mesh_list"])):
#my_current_mesh = my_input_cls['mesh_list'][0]
my_current_mesh = my_input_cls['mesh_list'][current_mesh]
#step-1 : Convert current_mesh_name into a char_list
#step-2: Select the 4th-to-last element in the char_list
#Step-3: Convert the 4th-to-last element to integer
#Now I can run any single mesh I want in whichever order
mesh_number=0 #For saving files according to mesh number
mesh_number = int(list(my_current_mesh)[-4])
print(mesh_number)
#MESH
my_mesh_cls = H5Mesh(my_current_mesh) #<-- Class instance is created
#print(my_current_mesh)
print("Maximum value of the mesh: ", my_mesh_cls.mesh.hmax())
#df.plot(my_mesh_cls.mesh)
#FUNCTION SPACE
problem_type = my_input_cls['problem_type']
number_of_moments = my_input_cls['number_of_moments']
my_function_space_cls = FunctionSpace(my_input_cls, my_mesh_cls) #<-- class instance
my_function_space_cls.set_function_space()
#SYSTEM MATRIX
my_system_matrix_cls = SystemMatrix(my_input_cls, my_mesh_cls) # <-- class instance
my_system_matrix_cls.convert_to_ufl_form()
#VARIATIONAL PROBLEM
my_var_prob_cls = FormulateVariationalProblem(
my_input_cls,
my_mesh_cls,
my_system_matrix_cls,
my_function_space_cls
) #<-- class instance
my_var_prob_cls.create_lhs()
my_var_prob_cls.create_rhs()
#SOLVER
my_solver_cls = Solver(my_input_cls, my_function_space_cls, my_var_prob_cls) #<-- class instance
if problem_type == 'nonlinear':
my_solver_cls.inbuilt_newton_solver()
u = my_solver_cls.u
else:
my_solver_cls.inbuilt_linear_solver()
u = my_solver_cls.u_Function
#===============
#Post-processing
#===============
sol = u.split()
def write_func(field_name, variable_name, mesh_num):
xdmffile_u = df.XDMFFile(df.MPI.comm_world,
'results_mathematica/{0}_{1}.xdmf'
.format(variable_name,mesh_num))
xdmffile_u.write(field_name)
xdmffile_u.close()
for i in range(my_input_cls["number_of_moments"]):
write_func(sol[i], i, mesh_number)
print("Program terminated successfully")
moment_order = my_input_cls["moment_order"]
def ErrorCalculation(exact_sol,numerical_sol):
#interpolate on the mesh
es = df.interpolate(exact_sol,my_function_space_cls.V.sub(0).collapse())
ns = df.interpolate(numerical_sol,my_function_space_cls.V.sub(0).collapse())
#compute values at the vertex
err_linf = np.max(np.abs(es.compute_vertex_values()- ns.compute_vertex_values()))
#err_l2 = np.linalg.norm(es.compute_vertex_values()-ns.compute_vertex_values())
err_l2 = df.errornorm(es,ns,"L2")#fenics inbuilt norm calculator
max_l2= np.linalg.norm(es.compute_vertex_values())
max_linf = np.max(np.abs(es.compute_vertex_values())) or 1
normalised_err_l2= err_l2/max_linf
normalised_err_linf = err_linf/max_linf
return normalised_err_l2, normalised_err_linf
def PressureErrorCalculation(exact_sol,numerical_sol1,numerical_sol2):
#interpolate on the mesh
es = interpolate(exact_sol, my_function_space_cls.V.sub(0).collapse())
ns1 = interpolate(numerical_sol1, my_function_space_cls.V.sub(0).collapse())
ns2 = interpolate(numerical_sol2, my_function_space_cls.V.sub(0).collapse())
#compute values at the vertex
#err_l2 = np.linalg.norm(es.compute_vertex_values()-
# (ns1.compute_vertex_values()+ns2.compute_vertex_values()))
err_l2 = df.errornorm(es,ns1+ns2,"L2")#fenics inbuilt norm calculator
err_linf = np.max(np.abs(es.compute_vertex_values()-
(ns1.compute_vertex_values()+ns2.compute_vertex_values())))
max_l2= np.linalg.norm(es.compute_vertex_values())
max_linf = np.max(np.abs(es.compute_vertex_values())) or 1
normalised_err_l2= err_l2/max_linf
normalised_err_linf = err_linf/max_linf
return normalised_err_l2, normalised_err_linf
if my_input_cls['moment_order']== 3:
with open("01_coeffs.cpp", "r") as file:
exact_solution_cpp_code = file.read()
load_value = df.compile_cpp_code(exact_solution_cpp_code)
#Temperature
t_e = df.CompiledExpression(load_value.Temperature(),degree=2)
t_l2,t_linf = ErrorCalculation(t_e,sol[0])
#Heat flux
sx_e = df.CompiledExpression(load_value.Heatfluxx(),degree=2)
sx_l2,sx_linf = ErrorCalculation(sx_e,sol[1])
sy_e = df.CompiledExpression(load_value.Heatfluxy(),degree=2)
sy_l2,sy_linf = ErrorCalculation(sy_e,sol[2])
errors =[
t_l2, t_linf,
sx_l2, sx_linf,
sy_l2,sy_linf ,
]
print(errors)
if my_input_cls['moment_order']== 'nono-6':
with open("01_coeffs.cpp", "r") as file:
exact_solution_cpp_code = file.read()
load_value = compile_cpp_code(exact_solution_cpp_code)
#Pressure
#p_e = CompiledExpression(load_value.Pressure(),degree=2)
#p_l2, p_linf = PressureErrorCalculation(p_e,sol[0],sol[3])
#velocity
ux_e = CompiledExpression(load_value.Velocityx(),degree=2)
ux_l2,ux_linf = ErrorCalculation(ux_e,sol[1])
uy_e = CompiledExpression(load_value.Velocityy(),degree=2)
uy_l2,uy_linf = ErrorCalculation(uy_e,sol[2])
#Pressure
t_e = CompiledExpression(load_value.Pressure(),degree=2)
t_l2,t_linf = ErrorCalculation(t_e,sol[0])
#Stress
sxx_e = CompiledExpression(load_value.Stressxx(),degree=2)
sxx_l2,sxx_linf = ErrorCalculation(sxx_e,sol[4])
sxy_e = CompiledExpression(load_value.Stressxy(),degree=2)
sxy_l2,sxy_linf = ErrorCalculation(sxy_e,sol[5])
syy_e = CompiledExpression(load_value.Stressyy(),degree=2)
syy_l2,syy_linf = ErrorCalculation(syy_e,sol[6])
#Heat flux
#sx_e = CompiledExpression(load_value.Heatfluxx(),degree=2)
#sx_l2,sx_linf = ErrorCalculation(sx_e,sol[7])
#sy_e = CompiledExpression(load_value.Heatfluxy(),degree=2)
#sy_l2,sy_linf = ErrorCalculation(sy_e,sol[8])
errors =[
t_l2, t_linf,
ux_l2, ux_linf,
uy_l2,uy_linf ,
sxx_l2,sxx_linf,
sxy_l2,sxy_linf,
syy_l2,syy_linf
]
print(errors)
#%%%%%%%%%%%%%%%%% Heat system Debug norm calculatation
###########################################################
'''
#interpolating the function values on function space
def interpolate_func(exact_sol,numerical_sol):
field_exact = df.interpolate(exact_sol, my_function_space_cls.V.sub(0).collapse())
field_numeric = df.interpolate(numerical_sol,my_function_space_cls.V.sub(0).collapse())
return field_exact, field_numeric
##
#Calculating L2 and Linf error
def error_calc(field_exact,field_numerical):
max_field_exact = np.max(np.abs(field_exact.compute_vertex_values())) or 1
err_L2 = df.errornorm(field_exact,field_numerical,"L2")
err_linf = np.max(np.abs(
field_exact.compute_vertex_values()
- field_numerical.compute_vertex_values()))
return max_field_exact,err_L2,err_linf
##OLD CODE BLOCK TO READ EXACT SOLUTION FROM CPP FILES
with open("01_coeffs.cpp", "r") as file:
exact_solution_cpp_code = file.read()
load_value = df.compile_cpp_code(exact_solution_cpp_code)
pressure_exact = df.CompiledExpression(load_value.Temperature(),degree=2)
#sx_exact = CompiledExpression(load_value.VelocityX(),degree=2)
#sy_exact = CompiledExpression(load_value.VelocityY(),degree=2)
#pressure
field_e_p,field_p = interpolate_func(pressure_exact,sol[0])
#write_func(field_e_p,'e_p')
max_exact_p,p_L2,p_linf = error_calc(field_e_p,field_p)
print("Theta L2,Linf :",p_L2/max_exact_p,p_linf/max_exact_p)
'''
examples/ffme/modules/FormulateVariationalProblem.py 0 → 100644
+ 263
− 0
View file @ 9693c9d3
import dolfin as df
import numpy
class FormulateVariationalProblem():
def __init__(self,
my_input_cls,
my_mesh_cls,
my_system_matrix_cls,
my_function_space_cls
):
#Input
self.my_input_cls = my_input_cls
self.my_system_matrix_cls = my_system_matrix_cls
self.my_mesh_cls = my_mesh_cls
self.my_function_space_cls = my_function_space_cls
#Class Input
self.moment_order = self.my_input_cls['moment_order']
self.number_of_moments = self.my_input_cls['number_of_moments']
self.problem_type = self.my_input_cls['problem_type']
self.Kn = self.my_input_cls['Kn']
self.Ma = self.my_input_cls['Ma']
self.stab_enable = self.my_input_cls['stabilization']['enable']
self.stab_type = self.my_input_cls['stabilization']['stab_type']
self.DELTA_T = self.my_input_cls['stabilization']['cip']['DELTA_T']
self.DELTA_P = self.my_input_cls['stabilization']['cip']['DELTA_P']
self.DELTA_U = self.my_input_cls['stabilization']['cip']['DELTA_U']
self.ht = self.my_input_cls['stabilization']['cip']['ht']
self.hp = self.my_input_cls['stabilization']['cip']['hp']
self.hu = self.my_input_cls['stabilization']['cip']['hu']
self.epsilon_w = self.my_input_cls['epsilon_w']
self.chi = self.my_input_cls['chi']
self.theta_w_inner = self.my_input_cls['bc'][3000]['theta_w']
self.u_t_w_inner = self.my_input_cls['bc'][3000]['u_t_w']
self.u_n_w_inner = self.my_input_cls['bc'][3000]['u_n_w']
self.p_w_inner = self.my_input_cls['bc'][3000]['p_w']
self.theta_w_outer = self.my_input_cls['bc'][3100]['theta_w']
self.u_t_w_outer = self.my_input_cls['bc'][3100]['u_t_w']
self.u_n_w_outer = self.my_input_cls['bc'][3100]['u_n_w']
self.p_w_outer = self.my_input_cls['bc'][3100]['p_w']
#Class SystemMatrix
self.SA_x = self.my_system_matrix_cls.SA_x_ufl
self.SA_y = self.my_system_matrix_cls.SA_y_ufl
self.SP_Coeff = self.my_system_matrix_cls.SP_Coeff_ufl
self.BC1 = self.my_system_matrix_cls.BC1_ufl
self.BC2 = self.my_system_matrix_cls.BC2_ufl
self.BC1_rhs = self.my_system_matrix_cls.BC1_rhs_ufl
self.BC2_rhs = self.my_system_matrix_cls.BC2_rhs_ufl
self.Symmetrizer= self.my_system_matrix_cls.Symmetrizer
#Class FunctionSpace
self.u = self.my_function_space_cls.u
self.v = self.my_function_space_cls.v
self.V = self.my_function_space_cls.V
#Class Mesh
self.mesh = self.my_mesh_cls.mesh
self.nv = self.my_mesh_cls.nv
self.ds = self.my_mesh_cls.ds
self.dS = self.my_mesh_cls.dS
#Output
self.a = 0.0
self.L = 0.0
self.F = 0.0
def create_a(self):
local_a =0.0
df.ds = self.ds #Import <----
local_a += (+0.5 * ((df.inner(self.v, (self.SA_x * df.Dx(self.u, 0)))) +\
(df.inner(self.v, (self.SA_y * df.Dx(self.u, 1))))) ) * df.dx
local_a += (-0.5 * ((df.inner(df.Dx(self.v, 0),(self.SA_x * self.u))) +\
(df.inner(df.Dx(self.v, 1),(self.SA_y * self.u)))) ) * df.dx
local_a += (+0.5 * df.inner(self.v, ((self.BC1 + self.BC2) * self.u)) ) * df.ds
local_a += (df.inner(self.v, (self.SP_Coeff * self.u)) ) * df.dx
return local_a
def apply_stabilization(self):
df.dS = self.dS #Important <----
h_msh = df.CellDiameter(self.mesh)
h_avg = ( h_msh("+") + h_msh("-") )/2.0
#Output
local_stab = 0.0
if self.stab_type == 'cip':
if self.moment_order==3 or self.moment_order=='ns':
local_stab += self.DELTA_T * h_avg**self.ht * df.jump(df.grad(self.u[0]),self.nv)\
* df.jump(df.grad(self.v[0]),self.nv) * df.dS #cip for temp
if self.moment_order==6:
local_stab += self.DELTA_P * h_avg**self.hp * df.jump(df.grad(self.u[0]),self.nv)\
* df.jump(df.grad(self.v[0]),self.nv) * df.dS #cip for pressure
local_stab += self.DELTA_U * h_avg**self.hu * df.jump(df.grad(self.u[1]),self.nv)\
* df.jump(df.grad(self.v[1]),self.nv) * df.dS #cip for velocity_x
local_stab += self.DELTA_U * h_avg**self.hu * df.jump(df.grad(self.u[2]),self.nv)\
* df.jump(df.grad(self.v[2]),self.nv) * df.dS #cip for velocity_y
if self.moment_order == 13:
local_stab += self.DELTA_T * h_avg**self.ht * df.jump(df.grad(self.u[3]),self.nv)\
* df.jump(df.grad(self.v[3]),self.nv) * df.dS #cip for temp
local_stab += self.DELTA_P * h_avg**self.hp * df.jump(df.grad(self.u[0]),self.nv)\
* df.jump(df.grad(self.v[0]),self.nv) * df.dS #cip for pressure
local_stab += self.DELTA_U * h_avg**self.hu * df.jump(df.grad(self.u[1]),self.nv)\
* df.jump(df.grad(self.v[1]),self.nv) * df.dS #cip for velocity_x
local_stab += self.DELTA_U * h_avg**self.hu * df.jump(df.grad(self.u[2]),self.nv)\
* df.jump(df.grad(self.v[2]),self.nv) * df.dS #cip for velocity_y
if self.moment_order == 'grad13':
local_stab += self.DELTA_T * h_avg**self.ht * df.jump(df.grad(self.u[3]),self.nv)\
* df.jump(df.grad(self.v[6]),self.nv) * df.dS #cip for temp
local_stab += self.DELTA_P * h_avg**self.hp * df.jump(df.grad(self.u[0]),self.nv)\
* df.jump(df.grad(self.v[0]),self.nv) * df.dS #cip for pressure
local_stab += self.DELTA_U * h_avg**self.hu * df.jump(df.grad(self.u[1]),self.nv)\
* df.jump(df.grad(self.v[1]),self.nv) * df.dS #cip for velocity_x
local_stab += self.DELTA_U * h_avg**self.hu * df.jump(df.grad(self.u[2]),self.nv)\
* df.jump(df.grad(self.v[2]),self.nv) * df.dS #cip for velocity_y
elif self.stab_type == 'gls':
local_stab = (0.0001* h_msh * df.inner(self.SA_x * df.Dx(self.u,0)\
+ self.SA_y * df.Dx(self.u,1) + self.SP_Coeff * self.u, self.SA_x * df.Dx(self.v,0)\
+ self.SA_y * df.Dx(self.v,1) + self.SP_Coeff * self.v ) * df.dx)
return local_stab
def inhomogeneity(self):
phi_local = df.Expression("atan2(x[1],x[0])",degree=2)
if self.moment_order == 3:
G_rhs_inner_list = [-2.0 * self.theta_w_inner, 0.0]
G_rhs_outer_list = [-2.0 * self.theta_w_outer, 0.0]
elif self.moment_order =='grad13':
#Input
u_n_w_inner = df.Expression("{}".format(self.u_n_w_inner),degree=2,phi=phi_local)
u_n_w_outer = df.Expression("{}".format(self.u_n_w_outer),degree=2,phi=phi_local)
u_t_w_inner = df.Expression("{}".format(self.u_t_w_inner),degree=2,phi=phi_local)
u_t_w_outer = df.Expression("{}".format(self.u_t_w_outer),degree=2,phi=phi_local)
#*******
G_rhs_inner_list = [-self.chi * u_n_w_inner,-self.chi * u_t_w_inner, 0.0]
G_rhs_outer_list = [-self.chi * u_n_w_outer,-self.chi * u_t_w_outer, 0.0]
elif self.moment_order == 'ns':
G_rhs_inner_list = [-1.0 *self.epsilon_w * self.theta_w_inner + self.u_n_w_inner\
,-1.0 * self.u_t_w_inner]
G_rhs_outer_list = [-1.0 *self.epsilon_w * self.theta_w_outer + self.u_n_w_outer\
,-1.0 * self.u_t_w_outer]
elif self.moment_order == 6:
#Input
u_n_w_inner = df.Expression("{}".format(self.u_n_w_inner),degree=2,phi=phi_local)
u_n_w_outer = df.Expression("{}".format(self.u_n_w_outer),degree=2,phi=phi_local)
u_t_w_inner = df.Expression("{}".format(self.u_t_w_inner),degree=2,phi=phi_local)
u_t_w_outer = df.Expression("{}".format(self.u_t_w_outer),degree=2,phi=phi_local)
p_w_inner = df.Expression("{}".format(self.p_w_inner),degree=2,phi=phi_local)
p_w_outer = df.Expression("{}".format(self.p_w_outer),degree=2,phi=phi_local)
#*******
G_rhs_inner_list = [-self.epsilon_w * self.chi * p_w_inner + u_n_w_inner,\
-self.chi * u_t_w_inner, 0.0, 0.0] # inner
G_rhs_outer_list = [-self.epsilon_w * self.chi * p_w_outer + u_n_w_outer,\
-self.chi * u_t_w_outer, 0.0, 0.0] # outer
elif self.moment_order == 13:
#Input
theta_w_inner = df.Expression("{}".format(self.theta_w_inner),degree=2,phi=phi_local)
theta_w_outer = df.Expression("{}".format(self.theta_w_outer),degree=2,phi=phi_local)
u_n_w_inner = df.Expression("{}".format(self.u_n_w_inner),degree=2,phi=phi_local)
u_n_w_outer = df.Expression("{}".format(self.u_n_w_outer),degree=2,phi=phi_local)
u_t_w_inner = df.Expression("{}".format(self.u_t_w_inner),degree=2,phi=phi_local)
u_t_w_outer = df.Expression("{}".format(self.u_t_w_outer),degree=2,phi=phi_local)
p_w_inner = df.Expression("{}".format(self.p_w_inner),degree=2,phi=phi_local)
p_w_outer = df.Expression("{}".format(self.p_w_outer),degree=2,phi=phi_local)
#*******
G_rhs_inner_list = [
- self.epsilon_w * self.chi * p_w_inner + u_n_w_inner,
- self.chi * u_t_w_inner,
- 2 * self.chi * theta_w_inner,
+ (2/5) * self.chi * theta_w_inner,
- (1/5) * self.chi * theta_w_inner,
+ self.chi * u_t_w_inner
] # inner
G_rhs_outer_list = [
- self.epsilon_w * self.chi * p_w_outer + u_n_w_outer,
- self.chi * u_t_w_outer,
- 2 * self.chi * theta_w_outer,
+ (2/5) * self.chi * theta_w_outer,
- (1/5) * self.chi * theta_w_outer,
+ self.chi * u_t_w_outer
] # outer
G_rhs_inner = df.as_vector(G_rhs_inner_list) #UFL form conversion
G_rhs_outer = df.as_vector(G_rhs_outer_list)
return G_rhs_inner, G_rhs_outer
def create_bc(self):
#Input
G_rhs_inner, G_rhs_outer = self.inhomogeneity()
local_bc = 0.0
df.ds = self.my_mesh_cls.ds #<---- Important
local_bc += (-0.5 * df.inner(self.v,((self.BC1_rhs-self.BC2_rhs) \
* G_rhs_inner)) ) * df.ds(3000) #RHS-2
local_bc += (-0.5 * df.inner(self.v,((self.BC1_rhs-self.BC2_rhs)\
* G_rhs_outer)) ) * df.ds(3100) #RHS-2
return local_bc
def add_nonlinearity(self):
#Output
local_NL =0.0
if self.moment_order=='ns':
local_NL = numpy.array(
[
df.Constant(0.0),
df.Constant(0.0),
df.Constant(0.0),
-((2.0/3.0) * self.u[1]**2 -(1.0/3.0) * self.u[2]**2),
-(self.u[1] * self.u[2]),
-(-(1.0/3.0) * self.u[1]**2 + (2.0/3.0) * self.u[2]**2)
]
)
local_NL = -(1/self.Kn)* local_NL
local_NL = numpy.dot(self.Symmetrizer,local_NL)
local_NL = df.as_vector(local_NL)
local_NL = df.inner(self.v,local_NL) * df.dx
elif self.moment_order == 'grad13':
local_NL = numpy.array(
[
df.Constant(0.0),
df.Constant(0.0),
df.Constant(0.0),
(self.Ma/self.Kn) * ((2.0/3.0) * self.u[1]**2 - (1.0/3.0) * self.u[2]**2),
(self.Ma/self.Kn) * (self.u[1] * self.u[2]),
(self.Ma/self.Kn) * (-(1.0/3.0) * self.u[1]**2 + (2.0/3.0) * self.u[2]**2),
df.Constant(0.0),
(self.Ma/self.Kn) * ( ( (10/9) * self.u[6] * self.u[1] )
+ ( (1/3) * (self.u[3]*self.u[1] + self.u[4] * self.u[2]) ) ),
(self.Ma/self.Kn) * ( ( (10/9) * self.u[6] * self.u[2] )
+ ( (1/3) * (self.u[4]*self.u[1] + self.u[5] * self.u[2]) ) )
]
)
local_NL = numpy.dot(self.Symmetrizer,local_NL)
local_NL = df.as_vector(local_NL)
local_NL = df.inner(self.v,local_NL) * df.dx
return local_NL
def source_term(self):
if self.moment_order==3:
F_rhs = [0.0]*self.number_of_moments
F_rhs[0] = df.Expression("2.0 - 1.0 * pow(sqrt(pow(x[0],2)+pow(x[1],2)),2)",degree=2)
elif self.moment_order== 'nono6':
F_rhs = [0.0]*NoV
elif self.moment_order== 'nono13':
F_rhs = [0.0]*NoV
R= Expression("sqrt(pow(x[0],2)+pow(x[1],2))",degree=2)
F_rhs[0] = Expression("1.0 * (1.0 - (5.0*pow(R,2))/(18.0*pow(kn,2))) * cos(phi)",R=R,kn=Kn,phi=phi,degree=2)
F_rhs[3] = Expression("1.0 * (1.0 - (5.0*pow(R,2))/(18.0*pow(kn,2))) * cos(phi)",R=R,kn=Kn,phi=phi,degree=2)
F_rhs[0] = Expression("2.0 - 1.0 * pow(sqrt(pow(x[0],2)+pow(x[1],2)),2)",degree=2)
F_rhs[3] = Expression("2.0 - 1.0 * pow(sqrt(pow(x[0],2)+pow(x[1],2)),2)",degree=2)
if self.moment_order == 3:
F_rhs = df.as_vector(F_rhs) # converting to ufl vector
local_source_term = df.inner(self.v, F_rhs) * df.dx
return local_source_term
#Left Hand Side
def create_lhs(self):
self.a = self.create_a()
if self.stab_enable == True:
if self.stab_type=='cip':
self.a += self.apply_stabilization()
elif self.stab_type == 'gls':
self.a += df.lhs(self.apply_stabilization())
def create_rhs(self):
self.L = self.create_bc()
if self.moment_order == 3:
self.L += self.source_term()
if self.stab_type == 'gls':
self.L += df.rhs(self.apply_stabilization())
if self.problem_type == 'nonlinear':
self.L += self.add_nonlinearity()
self.F = self.a - self.L
examples/ffme/modules/FunctionSpace.py 0 → 100644
+ 32
− 0
View file @ 9693c9d3
import dolfin as df
class FunctionSpace():
"""FunctionSpace."""
def __init__(self, my_input_cls, my_mesh_cls):
"""__init__.
:param my_input_cls:
:param my_mesh_cls:
"""
#Inputs
self.mesh = my_mesh_cls.mesh
self.number_of_moment = my_input_cls['number_of_moments']
self.problem_type = my_input_cls['problem_type']
#Returns
self.u = 0
self.v = 0
self.V = 0
def set_function_space(self):
"""set_function_space."""
self.V = df.VectorFunctionSpace(self.mesh,'P',1,dim=self.number_of_moment)
#Set test function(s)
v_list = df.TestFunctions(self.V)
#Convert to ufl form
self.v = df.as_vector(v_list)
#set trial function(s)
if self.problem_type == 'nonlinear':
u_list = df.Function(self.V)
elif self.problem_type == 'linear':
u_list = df.TrialFunctions(self.V)
#Convert to ufl form
self.u = df.as_vector(u_list)
examples/ffme/modules/Input.py 0 → 100644
+ 4
− 0
View file @ 9693c9d3
import yaml
class Input:
def __init__(self, input_file):
return 0
examples/ffme/modules/Mesh.py 0 → 100644
+ 23
− 0
View file @ 9693c9d3
import dolfin as df
class H5Mesh:
"""H5Mesh."""
def __init__(self, h5_file):
"""__init__.
:param h5_file:
"""
self.mesh = df.Mesh()
hdf = df.HDF5File(self.mesh.mpi_comm(),h5_file,"r")
hdf.read(self.mesh,"/mesh", False)
dim = self.mesh.topology().dim()
self.subdomains = df.MeshFunction("size_t", self.mesh, dim)
hdf.read(self.subdomains, "/subdomains")
self.boundaries = df.MeshFunction("size_t", self.mesh, dim-1)
hdf.read(self.boundaries, "/boundaries")
#Used in the program
self.nv = df.FacetNormal(self.mesh) #normal vector
self.ds = df.Measure("ds",domain=self.mesh, subdomain_data=self.boundaries)
self.dS = df.Measure("dS",domain=self.mesh, subdomain_data=self.boundaries)
examples/ffme/modules/Solver.py 0 → 100644
+ 71
− 0
View file @ 9693c9d3
import dolfin as df
class Solver():
def __init__(self,my_input_cls, my_function_space_cls, my_var_prob_cls):
self.my_input_cls = my_input_cls
self.my_var_prob_cls = my_var_prob_cls
self.my_function_space_cls = my_function_space_cls
self.u = my_function_space_cls.u
def inbuilt_newton_solver(self):
V = self.my_function_space_cls.V
F = self.my_var_prob_cls.F
self.u = self.my_var_prob_cls.u
du = df.TrialFunction(V) #TrialFunctions --> wrong, without 's'
Jac = df.derivative(F, self.u, du)
#user given solver parameters
abs_tol = self.my_input_cls['newton_abs_tol']
rel_tol = self.my_input_cls['newton_rel_tol']
max_itr = self.my_input_cls['newton_max_itr']
step_size = self.my_input_cls['newton_relaxation_parameter']
problem = df.NonlinearVariationalProblem(F,self.u,[],Jac)
solver = df.NonlinearVariationalSolver(problem)
solver.parameters ['newton_solver']['linear_solver'] = 'mumps'
solver.parameters ['newton_solver']['absolute_tolerance'] = abs_tol
solver.parameters ['newton_solver']['relative_tolerance'] = rel_tol
solver.parameters ['newton_solver']['maximum_iterations'] = max_itr
solver.parameters ['newton_solver']['relaxation_parameter'] =step_size
#df.solve(F == 0, self.u, [], J=Jac, solver_parameters={'newton_solver' : {'linear_solver' : 'mumps'}})
solver.solve()
#df.solve((a-L)==0, u,[])
def custom_newton_solver(self):
V = self.my_function_space_cls.V
F = self.my_var_prob_cls.F
#self.u = my_function_space_cls.u
du = df.TrialFunction(V) #TrialFunctions --> wrong, without 's'
Jac = df.derivative(F, self.u, du)
absolute_tol = 1E-5
relative_tol = 9E-2
u_inc = df.Function(V)
nIter = 0
relative_residual = 1
CONVERGED =True
MAX_ITER = 5000
while relative_residual > relative_tol and nIter < MAX_ITER and CONVERGED:
nIter += 1
A, b = df.assemble_system(Jac, -(F), [])
df.solve(A,u_inc.vector(),b)
#df.solve(F == 0, u_inc.vector(), [], J=Jac, solver_parameters={'newton_solver' : {'linear_solver' : 'mumps'}})
relative_residual = u_inc.vector().norm('l2')
#a = df.assemble(F)
residual = b.norm('l2')
lmbda = 0.8
self.u.vector()[:] += lmbda*u_inc.vector()
if residual > 1000:
CONVERGED = False
print("Did not converge after {} Iterations".format(nIter))
elif residual < absolute_tol:
CONVERGED =False
print("converged after {} Iterations".format(nIter))
print ('{0:2d} {1:3.2E} {2:5e}'.format(nIter, relative_residual, residual))
def inbuilt_linear_solver(self):
a = self.my_var_prob_cls.a
L = self.my_var_prob_cls.L
V = self.my_function_space_cls.V
self.u_Function = df.Function(V)
df.solve(a==L, self.u_Function,[],solver_parameters={'linear_solver':'mumps'})
examples/ffme/modules/SystemMatrix.py 0 → 100644
+ 75
− 0
View file @ 9693c9d3
#My library
from modules.matrices_dir import matrix_Ax
from modules.matrices_dir import matrix_Ay
from modules.matrices_dir import matrix_PCoeff
from modules.matrices_dir import matrix_Symm
from modules.matrices_dir import matrix_Tn
from modules.matrices_dir import matrix_Podd
from modules.matrices_dir import matrix_Peven
from modules.matrices_dir import matrix_BlockA
from modules.matrices_dir import matrix_L
import numpy as np
import dolfin as df
class SystemMatrix:
def __init__(self,my_input_cls, my_mesh_cls):
#Input
self.moment_order=my_input_cls['moment_order']
self.Kn = my_input_cls['Kn']
self.chi = my_input_cls['chi']
self.epsilon_w = my_input_cls['epsilon_w']
self.nv= my_mesh_cls.nv
def read_system(self):
A_x_list = (matrix_Ax.AxMatrix(self.moment_order))
A_y_list = (matrix_Ay.AyMatrix(self.moment_order))
P_Coeff_list = (matrix_PCoeff.PCoeffMatrix(self.moment_order,self.Kn))
Symmetrizer_list = (matrix_Symm.SymmMatrix(self.moment_order))
T_n_list = (matrix_Tn.TnMatrix(self.moment_order,self.nv))
p_odd_list = (matrix_Podd.PoddMatrix(self.moment_order))
p_even_list = (matrix_Peven.PevenMatrix(self.moment_order))
Block_A_list = (matrix_BlockA.BlockAMatrix(self.moment_order))
L_matrix_list = (matrix_L.LMatrix(self.moment_order,self.chi,self.epsilon_w))
self.A_x = np.array(A_x_list)
self.A_y = np.array(A_y_list)
self.P_Coeff = np.array(P_Coeff_list)
self.Symmetrizer = np.array(Symmetrizer_list)
self.T_n = np.array(T_n_list)
self.p_odd = np.array(p_odd_list)
self.p_even = np.array(p_even_list)
self.Block_A = np.array(Block_A_list)
self.L_matrix = np.array(L_matrix_list)
#Extra calculations
self.Block_A_trans = np.transpose(Block_A_list)
self.L_inverse = np.linalg.inv(L_matrix_list)
def compute_system(self):
self.SA_x = np.dot(self.Symmetrizer,self.A_x)
self.SA_y = np.dot(self.Symmetrizer,self.A_y)
self.SP_Coeff = np.dot(self.Symmetrizer,self.P_Coeff)
#Calculate matrices of boundary integral term
self.BC1 = np.dot(np.dot(np.dot(np.transpose(self.T_n),\
np.transpose(self.p_odd)),self.L_inverse),np.dot(self.p_odd,self.T_n))
self.BC2 = np.dot(np.dot(np.dot(np.dot(np.dot(np.transpose(self.T_n),\
np.transpose(self.p_even)),self.Block_A_trans),self.L_matrix),\
self.Block_A),np.dot(self.p_even,self.T_n))
self.BC1_rhs = np.dot(np.transpose(self.T_n),\
np.dot(np.transpose(self.p_even),self.Block_A_trans))
self.BC2_rhs = np.dot(np.transpose(self.T_n),\
np.dot(np.transpose(self.p_odd),self.L_inverse))
def convert_to_ufl_form(self):
#call the member functions
self.read_system()
self.compute_system()
#
self.SA_x_ufl = df.as_matrix(self.SA_x)
self.SA_y_ufl = df.as_matrix(self.SA_y)
self.SP_Coeff_ufl = df.as_matrix(self.SP_Coeff)
self.BC1_ufl = df.as_matrix(self.BC1)
self.BC2_ufl = df.as_matrix(self.BC2)
self.BC1_rhs_ufl = df.as_matrix(self.BC1_rhs)
self.BC2_rhs_ufl = df.as_matrix(self.BC2_rhs)
examples/ffme/modules/matrices_dir/matrix_Ax.py 0 → 100644
+ 95
− 0
View file @ 9693c9d3
#Store Ax matrices
def AxMatrix(moment_order):
if moment_order==3:
Ax = [
[0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
[5.0/2.0, 0.0, 0.0, 1.0/2.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 1.0/2.0, 0.0],
[0.0, 16.0/5.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 12.0/5.0, 0.0, 0.0, 0.0],
[0.0, -8.0/5.0, 0.0, 0.0, 0.0, 0.0]
]
return Ax
elif moment_order=='grad13':
Ax = [
[0, 1, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 1, 0, 0, 2/3, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0, 0],
[0, 4/3, 0, 0, 0, 0, 0, 8/15, 0],
[0, 0, 1, 0, 0, 0, 0, 0, 2/5],
[0, -(2/3), 0, 0, 0, 0, 0, -(4/15), 0],
[0, 1, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 1, 0, 0, 5/3, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0, 0]
]
return Ax
elif moment_order==6:
Ax = [
[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 4.0/3.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
[0.0, -2.0/3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
[0.0, 0.0, 0.0, 6.0/5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 16.0/15.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, -4.0/15.0, 0.0, 2.0/3.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, -4.0/5.0, 0.0, 0.0, 0.0, 0.0, 0.0]
]
return Ax
elif moment_order==13:
Ax = [
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 4/3, 0, 0, 0, 0, 0, 8/15, 0, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0, 2/5, 0, 1, 0, 0, 0, 0, 0, 0],
[0, -(2/3), 0, 0, 0, 0, 0, -(4/15), 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 5/2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1/6, 1/2, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0,0, 0, 0, 0, 0, 0, 1/2, 0],
[0, 0, 0, 0, 9/5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 8/5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -(2/5), 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0,0, 0, 0, -(6/5), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 56/15, 0,0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 14/5, 0, 0, 0, 0,0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, -(28/15), 0, 0, 0, 0, 0, 0, 0, 0, 0]
]
return Ax
elif moment_order=='ns':
Ax = [
[0, 1, 0, 0, 0, 0],
[1, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 1, 0],
[0, 2/3,0, 0, 0, 0],
[0, 0, 1/2, 0, 0, 0],
[0, -(1/3), 0, 0, 0, 0]
]
return Ax
elif moment_order=='g26w':
Ax = [
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, -0.6666666666666666, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1.3333333333333333, 0, 0, 0, 0, 0, -0.5333333333333333, 0, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0, -0.4, 0, 1, 0, 0, 0, 0, 0, 0],
[0, -0.6666666666666666, 0, 0, 0, 0, 0, 0.26666666666666666, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1.6666666666666667, -1, 0, 0, 0, 0, 0, 0, 0, 0, -0.6666666666666666, 1, 0, 0],
[0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 1.8, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5142857142857142, 0, 0],
[0, 0, 0, 0, 0, 1.6, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.45714285714285713, 0],
[0, 0, 0, 0, -0.4, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.11428571428571428, 0, -0.2857142857142857],
[0, 0, 0, 0, 0, -1.2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.34285714285714286, 0],
[0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1.8666666666666667, 0, -1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 1.4, 0, -1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, -0.9333333333333333, 0, 0, 0, -1, 0, 0, 0, 0, 0]
]
return Ax
else :
print("Wrong moment number")
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment