Skip to content
Snippets Groups Projects
Commit 9cccc616 authored by David Seus's avatar David Seus
Browse files

restructure repo

parent f254391c
No related branches found
No related tags found
No related merge requests found
Showing with 275 additions and 115 deletions
#!/usr/bin/python3
# I have been struggling about having a list of node indices:
#
# nodes = [1, 10, 22, ..., 32]
#
# of nodes at the boundary and an associated list of values
#
# node_values = [10.2, 1.2,0.1, ..., -0.2]
#
# and how to assign this to a function defined on the mesh, in order to use it as boundary conditions, instead of using an analytical function or a constant.
#
# After searchin thoroughly I manage to find a solution that works for low order basis functions, for both scalars and vector. I put here my solution to see if this is the way to do it and how can I extend it to higher order basis functions. Also, if it helps someone:
import numpy
from time import time
from dolfin import *
def dofs_to_verts(mesh,boundaryFunction,boundaryValue):
"""
Compute the mapping from vertices to degrees of freedom only at the boundary.
This saves time, since we do not have to loop over all the cells of the domain.
If the for loop was replaced by a simple c++ function call, it would be much faster.
INPUTS ::
mesh :: The mesh where we want to compute the mapping between vertices
and degrees of freedom.
boundaryFunction :: A MeshFunction which assigns integer values to
edges of the mesh.
boundaryValue :: An integer value which we use to identify the boundary
edges we want to use. The MeshFunction should have
this value assigned to some of the edges.
"""
# get the connectivity of the mesh
mesh.init(1) # initialize edges
mesh_topology = mesh.topology() # get the mesh topology
conn12 = mesh_topology(1,2) # get the cells to which an edge belongs
# get the number of edges
n_edges = boundaryFunction.size()
# get the indices of the boundary edges, as given by the boundaryFunction
boundaryEdges = numpy.arange(n_edges)
boundaryEdges = boundaryEdges[boundaryFunction.array() == boundaryValue]
# Define the function spaces associated to the mesh
# for now only CG of order 1 works for sure
V = FunctionSpace(mesh, "CG", 1)
Vv = VectorFunctionSpace(mesh, "CG", 1, dim=2)
# Allocate memory space of the array that stores the mapping from
# vertex number to degree of freedom number for both scalar and vector
# valued functions (note that for vector valued functions we need to indices
# since the x component and the y component have different indices)
vert_to_dofs = numpy.zeros(mesh.num_vertices(), dtype=numpy.uintp)
vectordofs_to_vert = numpy.zeros([2,mesh.num_vertices()], dtype=numpy.uintp)
# Get the degree of freedom maps
dm = V.dofmap()
dms = [Vv.sub(i).dofmap() for i in range(2)] # the i=0 refers to x component
# the i=1 refers to y component
# get the cells of the mesh
mesh_cells = mesh.cells()
# Since we are only interested in the boundary vertices we just loop over
# the boundary edges as given by the MeshFunction and stored in boundaryEdges
# each time we loop we get the numbering of the degrees of freedom of the
# vertices of the cell that contains the boundary edge and we construct
# the mapping from vertex number to degree of freedom number
for edge in boundaryEdges:
# get the cell index to which the edge belongs
cell_ind = conn12(edge)[0]
# get the vertices of the cell
vert_inds = mesh_cells[cell_ind]
# add the degree of freedom numbering to the vertices
vert_to_dofs[vert_inds] = numpy.array(dm.cell_dofs(cell_ind),dtype=numpy.uintp)
# now we do the same but for vector valued quantites
# since, in this case, our vector have two quantities (the x and y components)
# we have to iterate over the two sub dof maps
for i, (dms_i, dmcs_i) in enumerate(zip(dms, dms)):
vectordofs_to_vert[i,vert_inds] = numpy.array(dms_i.cell_dofs(cell_ind),dtype=numpy.uintp)
# Return vertex to dof mapping
return vert_to_dofs, vectordofs_to_vert
# this is an example of how to assign values at the nodes for both a scalar
# and a vector valued function to a function
# create a unit square mesh (any other mesh works)
mesh = UnitSquareMesh(20,20)
# --------------------------------------------------------------------------
# this part is optional, if you already know the indices of the vertices at
# the boundary just skip this part
# define the boundary mesh
boundaryMesh = BoundaryMesh(mesh)
# if you already don't have the boundary nodes indices
# get the boundary nodes indices
nodes = boundaryMesh.vertex_map().array()
# --------------------------------------------------------------------------
# --------------------------------------------------------------------------
# This part is just to generate a boundary function which assigns a value
# boundary_value to the boundary edges
boundary_value = 10
# get the global number of the boundary edges
edge_map = boundaryMesh.cell_map().array()
# get the number of edges at the boundary
n_boundaryEdges = edge_map.size
# define the boundary function which defines the boundary
boundaryFunction = MeshFunction('int',mesh,1)
n_edges = boundaryFunction.array().size # the total number of edges
# allocate memory space to boundary values
boundary_values = numpy.zeros(n_edges,dtype='uintp')
# assign the boundary_value to the array of boundary_values
boundary_values[edge_map] = boundary_value
# assign the boundary_values to the boundary_function
boundaryFunction.set_values(boundary_values)
# --------------------------------------------------------------------------
# generate some random values for the values to assign to the scalar function at the nodes
node_values_scalar = numpy.random.rand(nodes.size)
# generate some random values for the values to assign to the vector function at the nodes
node_values_vector = numpy.random.rand(2,nodes.size)
start = time()
# get the mapping of vertices at boundary to degrees of freedom
vert_to_dofs, vectordofs_to_vert = dofs_to_verts(mesh,boundaryFunction,boundary_value)
#print("It took :: %f s" %(time()-start)
# assign vertex values at the boundary for the scalar function
# generate the function space
V = FunctionSpace(mesh,'Lagrange',1)
# generate the scalar function
myFunction = Function(V)
# get the vector containing the data values at the degrees of freedom
myFunction_vector = myFunction.vector()
# use the mapping from vertices (nodes) to degrees of freedom to assign the
# vertex values of the scalar function
myFunction_vector[vert_to_dofs[nodes]] = node_values_scalar
plot(myFunction)
# assign vertex values at the boundary for the vector function
# generate the function space
Vv = VectorFunctionSpace(mesh, "Lagrange", 1, dim=2)
# generate the vector function
myVectorFunction = Function(Vv)
# get the vector containing the data values at the degrees of freedom
myVectorFunction_vector = myVectorFunction.vector()
# use the mapping from vertices (nodes) to degrees of freedom to assign the
# vertex values of the scalar function
myVectorFunction_vector[vectordofs_to_vert[0,nodes]] = node_values_vector[0,:] # x component
myVectorFunction_vector[vectordofs_to_vert[1,nodes]] = node_values_vector[1,:] # y component
plot(myVectorFunction)
interactive()
File moved
File moved
File moved
File moved
# Meshes
The mesh data in the various folders are only saved here to be able to plot the
data either with `paraview` or the script found in the folder `Plotscript`.
They are produced by the simulation class as output and are not reused.
Note that the scripts in `Plotscript` are courtesy of Jim Magiera.
# Numerical experiments for LDD-TPR
This project provides a Fenics based 2D domain decomposition code implementing an LDD solver for two-phase flow systems.
Two-phase flow systems means, that each subdomain can either adopt the full two-phase flow equations or the Richards/Richardson.
Flexible domain substructurings are possible so long as each subdomain has polygonal boundary.
## Manual how to set up latest fenix image in Docker
Pull the latest fenics Docker image
docker pull quay.io/fenicsproject/stable:latest
## Create docker container named LDD-TPR with graphical display, share folder `pwd`
docker run --dns=129.69.252.252 -it --env HOST_UID=$(id -u $USER) --env HOST_GID=$(id -g $USER) --env DISPLAY=unix$DISPLAY --device /dev/dri --volume /tmp/.X11-unix:/tmp/.X11-unix:rw --volume $(pwd):/home/fenics/shared --name LDD-TPR-fenics2019 quay.io/fenicsproject/stable:latest /bin/bash
If you want the container do be automatically deleted after exiting add `--rm` option.
If you have trouble with internet connection inside the container use
--dns=YOUR_DNS_SERVER
## Start LDD-TPR container and step into the container
docker start LDD-TPR & docker exec -ti -u fenics LDD-TPR /bin/bash -l
## Usefull docker commands
List all docker container
docker ps -a
Remove container
docker rm $container
Remove all stopped container
docker container prune
List all images
docker images
Remove image
docker rmi $image
# Troubleshooting
**Problem**
1. I can't install packages via `apt-get install` inside the container
2. I cannot create files or folders (no write permissions)
**Solution**
1. If the package is not found first `apt-get update` and then try again. If there is no connection
check your dns settings.
2. In the container, execute once the script `Rechtesetup/setpermissions.sh` to gain write access in `/home/fenics/shared`
cd /home/fenics/shared/Rechtesetup & sudo ./setpermissions.sh
## Nützliche FEniCS Links
- Forum, viele Fragen wurden evtl. schon gestellt, bzw. können ggf. gestellt werden:
https://www.allanswered.com/community/s/fenics-project/
- Python API reference (die von der neuesten Version existiert quasi noch nicht):
https://fenicsproject.org/docs/dolfin/2017.2.0/python/programmers-reference/index.html
- Fenics git wenn man wirklich mal in die source files schauen möchte:
https://bitbucket.org/fenics-project/
File moved
File moved
# Helperscripts
The script
merge_spacetime_errornorms.py
is part of the `../LDDsimulation/helpers.py` but has been adapted to manually
use it to merge the spacetime errornorm files that come out of mesh studies,
manually.
# global domain
subdomain0_vertices = [df.Point(0.0,0.0), #
df.Point(13.0,0.0),#
df.Point(13.0,8.0),#
df.Point(0.0,8.0)]
interface12_vertices = [df.Point(0.0, 7.0),
df.Point(9.0, 7.0),
df.Point(10.5, 7.5),
df.Point(12.0, 7.0),
df.Point(13.0, 6.5)]
# subdomain1.
subdomain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
interface12_vertices[2],
interface12_vertices[3],
interface12_vertices[4], # southern boundary, 12 interface
subdomain0_vertices[2], # eastern boundary, outer boundary
subdomain0_vertices[3]] # northern boundary, outer on_boundary
# vertex coordinates of the outer boundaries. If it can not be specified as a
# polygon, use an entry per boundary polygon. This information is used for defining
# the Dirichlet boundary conditions. If a domain is completely internal, the
# dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
0: [interface12_vertices[4], #
subdomain0_vertices[2], # eastern boundary, outer boundary
subdomain0_vertices[3],
interface12_vertices[0]]
}
# interface23
interface23_vertices = [df.Point(0.0, 5.0),
df.Point(3.0, 5.0),
# df.Point(6.5, 4.5),
df.Point(6.5, 5.0),
df.Point(9.5, 5.0),
# df.Point(11.5, 3.5),
# df.Point(13.0, 3)
df.Point(11.5, 5.0),
df.Point(13.0, 5.0)
]
#subdomain1
subdomain2_vertices = [interface23_vertices[0],
interface23_vertices[1],
interface23_vertices[2],
interface23_vertices[3],
interface23_vertices[4],
interface23_vertices[5], # southern boundary, 23 interface
subdomain1_vertices[4], # eastern boundary, outer boundary
subdomain1_vertices[3],
subdomain1_vertices[2],
subdomain1_vertices[1],
subdomain1_vertices[0] ] # northern boundary, 12 interface
subdomain2_outer_boundary_verts = {
0: [interface23_vertices[5],
subdomain1_vertices[4]],
1: [subdomain1_vertices[0],
interface23_vertices[0]]
}
# interface34
interface34_vertices = [df.Point(0.0, 2.0),
df.Point(4.0, 2.0),
df.Point(9.0, 2.5),
df.Point(10.5, 2.0),
df.Point(13.0, 1.5)]
# subdomain3
subdomain3_vertices = [interface34_vertices[0],
interface34_vertices[1],
interface34_vertices[2],
interface34_vertices[3],
interface34_vertices[4], # southern boundary, 34 interface
subdomain2_vertices[5], # eastern boundary, outer boundary
subdomain2_vertices[4],
subdomain2_vertices[3],
subdomain2_vertices[2],
subdomain2_vertices[1],
subdomain2_vertices[0] ] # northern boundary, 23 interface
subdomain3_outer_boundary_verts = {
0: [interface34_vertices[4],
subdomain2_vertices[5]],
1: [subdomain2_vertices[0],
interface34_vertices[0]]
}
# subdomain4
subdomain4_vertices = [subdomain0_vertices[0],
subdomain0_vertices[1], # southern boundary, outer boundary
subdomain3_vertices[4],# eastern boundary, outer boundary
subdomain3_vertices[3],
subdomain3_vertices[2],
subdomain3_vertices[1],
subdomain3_vertices[0] ] # northern boundary, 34 interface
subdomain4_outer_boundary_verts = {
0: [subdomain4_vertices[6],
subdomain4_vertices[0],
subdomain4_vertices[1],
subdomain4_vertices[2]]
}
subdomain_def_points = [subdomain0_vertices,#
subdomain1_vertices,#
subdomain2_vertices,#
subdomain3_vertices,#
subdomain4_vertices
]
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment