Skip to content

Commit

Permalink
housekeeping
Browse files Browse the repository at this point in the history
  • Loading branch information
pierricmora committed Jun 8, 2023
1 parent 24d4390 commit 17c566a
Show file tree
Hide file tree
Showing 22 changed files with 345 additions and 177 deletions.
8 changes: 7 additions & 1 deletion elastodynamicsx/__init__.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
# Copyright (C) 2023 Pierric Mora
#
# This file is part of ElastodynamiCSx
#
# SPDX-License-Identifier: MIT

"""Main module for ElastodynamiCSx"""

__version__ = '0.15'
__version__ = '0.2'
9 changes: 6 additions & 3 deletions elastodynamicsx/pde/__init__.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
# Copyright (C) 2023 Pierric Mora
#
# This file is part of ElastodynamiCSx
#
# SPDX-License-Identifier: MIT

"""
Tools for building a Partial Differential Equation from material laws. The
package also provides tools for building Boundary Conditions.
Expand All @@ -24,12 +30,9 @@
pde = PDE(materials=[mat1, mat2, ...], bodyforces=[bf1, bf2, ...])
"""


from .boundarycondition import *

from .pde import *
from .materials import *
#from .elasticmaterial import *
#from .hyperelasticmaterial import *
from .bodyforce import *
from .timeschemes import *
14 changes: 10 additions & 4 deletions elastodynamicsx/pde/bodyforce.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
# Copyright (C) 2023 Pierric Mora
#
# This file is part of ElastodynamiCSx
#
# SPDX-License-Identifier: MIT

import ufl

from elastodynamicsx.utils import get_functionspace_tags_marker
Expand All @@ -9,20 +15,20 @@ class BodyForce():
Representation of the rhs term (the 'F' term) of a pde such as defined
in the PDE class. An instance represents a single source.
"""

def __init__(self, functionspace_tags_marker, value, **kwargs):
self._value = value
function_space, cell_tags, marker = get_functionspace_tags_marker(functionspace_tags_marker)
md = kwargs.get('metadata', PDE.default_metadata)
self._dx = ufl.Measure("dx", domain=function_space.mesh, subdomain_data=cell_tags, metadata=md)(marker)



@property
def L(self):
"""The linear form function"""
return lambda v: ufl.inner(self._value, v) * self._dx


@property
def value(self):
return self._value


99 changes: 57 additions & 42 deletions elastodynamicsx/pde/boundarycondition.py
Original file line number Diff line number Diff line change
@@ -1,25 +1,32 @@
# Copyright (C) 2023 Pierric Mora
#
# This file is part of ElastodynamiCSx
#
# SPDX-License-Identifier: MIT

import numpy as np
from dolfinx import fem
from petsc4py import PETSc
import ufl

from elastodynamicsx.utils import get_functionspace_tags_marker


class BoundaryCondition():
"""
Representation of a variety of boundary conditions
Examples of use:
##### #####
# Free or Clamp #
##### #####
# Imposes sigma(u).n=0 on boundary n°1
bc = BoundaryCondition((V, facet_tags, 1), 'Free' )
# Imposes u=0 on boundaries n°1 and 2
bc = BoundaryCondition((V, facet_tags, (1,2)), 'Clamp')
# Apply BC to all boundaries when facet_tags and marker are not specified
# or set to None
bc = BoundaryCondition(V, 'Clamp')
Expand All @@ -28,21 +35,21 @@ class BoundaryCondition():
##### #####
# General #
##### #####
# Imposes u=u_D on boundary n°1
bc = BoundaryCondition((V, facet_tags, 1), 'Dirichlet', u_D)
# Imposes sigma(u).n=T_N on boundary n°1
bc = BoundaryCondition((V, facet_tags, 1), 'Neumann', T_N)
# Imposes sigma(u).n=r*(u-s) on boundary n°1
bc = BoundaryCondition((V, facet_tags, 1), 'Robin', (r, s))
##### #####
# Periodic #
##### #####
# Given: x(2) = x(1) - P
# Imposes: u(2) = u(1)
# Where: x(i) are the coordinates on boundaries n°1,2
Expand All @@ -55,30 +62,30 @@ class BoundaryCondition():
Px, Py, Pz = 0, 0, depth #for z-periodic, and tag=back
P = [Px,Py,Pz]
bc = BoundaryCondition((V, facet_tags, 2), 'Periodic', P)
##### #####
# BCs involving the velocity #
##### #####
# (for a scalar function_space)
# Imposes sigma(u).n = z*v on boundary n°1, with v=du/dt. Usually z=rho*c
bc = BoundaryCondition((V, facet_tags, 1), 'Dashpot', z)
# (for a vector function_space)
# Imposes sigma(u).n = z_N*v_N+z_T*v_T on boundary n°1,
# with v=du/dt, v_N=(v.n)n, v_T=v-v_N
# Usually z_N=rho*c_L and z_T=rho*c_S
bc = BoundaryCondition((V, facet_tags, 1), 'Dashpot', (z_N, z_T))
Adapted from:
https://jsdokken.com/dolfinx-tutorial/chapter3/robin_neumann_dirichlet.html
"""

### ### ### ###
### static ###
### ### ### ###

def get_dirichlet_BCs(bcs):
out = []
for bc in bcs:
Expand All @@ -87,14 +94,16 @@ def get_dirichlet_BCs(bcs):
if type(bc) == BoundaryCondition and issubclass(type(bc.bc), fem.DirichletBCMetaClass):
out.append(bc.bc)
return out



def get_mpc_BCs(bcs):
out = []
for bc in bcs:
if bc.type == 'periodic':
out.append(bc)
return out



def get_weak_BCs(bcs):
out = []
for bc in bcs:
Expand All @@ -104,62 +113,66 @@ def get_weak_BCs(bcs):
if not(test_mpc or test_d):
out.append(bc)
return out





### ### ### ### ###
### non-static ###
### ### ### ### ###



def __init__(self, functionspace_tags_marker, type_, values=None, **kwargs):
from elastodynamicsx.pde import PDE
#
function_space, facet_tags, marker = get_functionspace_tags_marker(functionspace_tags_marker)

type_ = type_.lower()
nbcomps = function_space.element.num_sub_elements #number of components if vector space, 0 if scalar space
nbcomps = function_space.element.num_sub_elements # number of components if vector space, 0 if scalar space

# shortcuts: Free->Neumann with value=0; Clamp->Dirichlet with value=0
if type_ == "free":
type_ = "neumann"
z0s = np.array([0]*nbcomps, dtype=PETSc.ScalarType) if nbcomps>0 else PETSc.ScalarType(0)
values = fem.Constant(function_space.mesh, z0s)

elif type_ == "clamp":
type_ = "dirichlet"
z0s = np.array([0]*nbcomps, dtype=PETSc.ScalarType) if nbcomps>0 else PETSc.ScalarType(0)
values = fem.Constant(function_space.mesh, z0s)

self._type = type_
self._values = values
md = kwargs.get('metadata', PDE.default_metadata)
ds = ufl.Measure("ds", domain=function_space.mesh, subdomain_data=facet_tags, metadata=md)(marker) #also valid if facet_tags or marker are None
ds = ufl.Measure("ds",
domain=function_space.mesh,
subdomain_data=facet_tags,
metadata=md)(marker) # also valid if facet_tags or marker are None

if type_ == "dirichlet":
fdim = function_space.mesh.topology.dim - 1
#facets = facet_tags.find(marker) #not available on dolfinx v0.4.1
#facets = facet_tags.find(marker) # not available on dolfinx v0.4.1
facets = facet_tags.indices[ facet_tags.values == marker ]
dofs = fem.locate_dofs_topological(function_space, fdim, facets)
self._bc = fem.dirichletbc(values, dofs, function_space) #fem.dirichletbc
self._bc = fem.dirichletbc(values, dofs, function_space) # fem.dirichletbc

elif type_ == "neumann":
self._bc = lambda v : ufl.inner(values, v) * ds #Linear form
self._bc = lambda v : ufl.inner(values, v) * ds # Linear form

elif type_ == "robin":
r_, s_ = values
if nbcomps>0:
r_, s_ = ufl.as_matrix(r_), ufl.as_vector(s_)
self._bc = lambda u,v: ufl.inner(r_*(u-s_), v)* ds #Bilinear form
self._bc = lambda u,v: ufl.inner(r_*(u-s_), v)* ds # Bilinear form

elif type_ == 'dashpot':
if nbcomps == 0: #scalar function space
self._bc = lambda u_t,v: values * ufl.inner(u_t, v)* ds #Bilinear form, to be applied on du/dt
else: #vector function space
if nbcomps == 0: # scalar function space
self._bc = lambda u_t,v: values * ufl.inner(u_t, v)* ds # Bilinear form, to be applied on du/dt
else: # vector function space
if function_space.mesh.topology.dim == 1:
self._bc = lambda u_t,v: ((values[0]-values[1])*ufl.inner(u_t[0], v[0]) + values[1]*ufl.inner(u_t, v))* ds #Bilinear form, to be applied on du/dt
self._bc = lambda u_t,v: ((values[0]-values[1])*ufl.inner(u_t[0], v[0]) + values[1]*ufl.inner(u_t, v))* ds # Bilinear form, to be applied on du/dt
else:
n = ufl.FacetNormal(function_space)
self._bc = lambda u_t,v: ((values[0]-values[1])*ufl.dot(u_t, n)*ufl.inner(n, v) + values[1]*ufl.inner(u_t, v))* ds #Bilinear form, to be applied on du/dt
self._bc = lambda u_t,v: ((values[0]-values[1])*ufl.dot(u_t, n)*ufl.inner(n, v) + values[1]*ufl.inner(u_t, v))* ds # Bilinear form, to be applied on du/dt

elif type_ == 'periodic-do-not-use': #TODO: try to calculate P from two given boundary markers
assert len(marker)==2, "Periodic BC requires two facet tags"
fdim = function_space.mesh.topology.dim - 1
Expand All @@ -180,20 +193,22 @@ def __init__(self, functionspace_tags_marker, type_, values=None, **kwargs):
def slave_to_master_map(x):
return x + P[:,np.newaxis]
self._bc = (facet_tags, marker, slave_to_master_map)

else:
raise TypeError("Unknown boundary condition: {0:s}".format(type_))


@property
def bc(self):
"""The boundary condition"""
return self._bc


@property
def type(self):
return self._type


@property
def values(self):
return self._values

6 changes: 6 additions & 0 deletions elastodynamicsx/pde/materials/__init__.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
# Copyright (C) 2023 Pierric Mora
#
# This file is part of ElastodynamiCSx
#
# SPDX-License-Identifier: MIT

from .material import *
from .elasticmaterial import *
from .anisotropicmaterials import *
Expand Down
Loading

0 comments on commit 17c566a

Please sign in to comment.