From 8c1fc8d097b6ce4b191dcf2aef02bf1c71f3da72 Mon Sep 17 00:00:00 2001 From: Gustafsson Tom Date: Tue, 10 Apr 2018 16:27:46 +0300 Subject: [PATCH] 3d elasticity example --- examples/ex11.py | 25 +++++++++++++++ examples/ex2.py | 2 +- examples/ex3.py | 2 +- examples/ex5.py | 2 +- examples/ex6.py | 2 +- setup.py | 2 +- skfem/assembly.py | 30 ++++++------------ skfem/mesh.py | 4 +-- skfem/models/elasticity.py | 65 ++++++++++++++++++++++++++++++++++++++ skfem/utils.py | 13 +++++++- 10 files changed, 119 insertions(+), 28 deletions(-) create mode 100644 examples/ex11.py diff --git a/examples/ex11.py b/examples/ex11.py new file mode 100644 index 000000000..53521d2aa --- /dev/null +++ b/examples/ex11.py @@ -0,0 +1,25 @@ +from skfem import * +from skfem.models.elasticity import * + +m = MeshHex() +m.refine(3) +ib = basis_elasticity(m) + +K = asm(linear_elasticity(*lame_parameters(1e3, 0.3)), ib) + +@boundary_clamped +def left(x, y, z): + return x==0.0 + +@boundary_prescribed(0.3, 0.0, 0.0) +def right(x, y, z): + return x==1.0 + +u, I = initialize(ib, left, right) + +u[I] = solve(*condense(K, 0*u, I=I, x=u)) + +sf = 1.0 +for itr in range(3): + m.p[itr, :] += sf*u[ib.dofnum.n_dof[itr, :]] +m.export_vtk('elasticity') diff --git a/examples/ex2.py b/examples/ex2.py index 128c14d37..c3127eda0 100644 --- a/examples/ex2.py +++ b/examples/ex2.py @@ -39,7 +39,7 @@ def ddot(T1,T2): K = asm(bilinf, ib) f = asm(unit_load, ib) -x, D = ib.essential_bc() +x, D = ib.find_dofs() I = ib.dofnum.complement_dofs(D) x[I] = solve(*condense(K, f, I=I), solver=solver_direct_cholmod()) diff --git a/examples/ex3.py b/examples/ex3.py index 88bf364f5..f44abbc2a 100644 --- a/examples/ex3.py +++ b/examples/ex3.py @@ -27,7 +27,7 @@ def mass(u, du, v, dv, w): M = asm(mass, gb) -y, D = gb.essential_bc(lambda x, y: x==0.0) +y, D = gb.find_dofs(lambda x, y: x==0.0) I = gb.dofnum.complement_dofs(D) diff --git a/examples/ex5.py b/examples/ex5.py index 32da69d87..b70764a8d 100644 --- a/examples/ex5.py +++ b/examples/ex5.py @@ -34,7 +34,7 @@ def facetlinf(v, dv, w): b = asm(facetlinf, fb) -_, D = ib.essential_bc(lambda x, y: (y == 0.0)) +_, D = ib.find_dofs(lambda x, y: (y == 0.0)) I = ib.dofnum.complement_dofs(D) import scipy.sparse diff --git a/examples/ex6.py b/examples/ex6.py index 55525d3b8..7ba6da01f 100644 --- a/examples/ex6.py +++ b/examples/ex6.py @@ -21,7 +21,7 @@ def linf(v, dv, w): f = asm(linf, ib) -x, D = ib.essential_bc() +x, D = ib.find_dofs() I = ib.dofnum.complement_dofs(D) x[I] = solve(*condense(K, f, D=D)) diff --git a/setup.py b/setup.py index 652f5d916..63df17e80 100644 --- a/setup.py +++ b/setup.py @@ -6,7 +6,7 @@ setup( name='scikit-fem', - version='0.1.3', + version='0.1.4', description='Simple finite element assemblers', long_description='Easy to use finite element assemblers and related tools. See Github page for more information and examples.', # Optional url='https://github.com/kinnala/scikit-fem', diff --git a/skfem/assembly.py b/skfem/assembly.py index a778eb0e0..38e072bb0 100644 --- a/skfem/assembly.py +++ b/skfem/assembly.py @@ -97,14 +97,12 @@ def interpolate(self, w, derivative=False): return W, dW return W - def essential_bc(self, test=None, bc=None, boundary=True, dofrows=None, check_vertices=True, - check_facets=True, check_edges=True): - """Helper function for setting essential boundary conditions. + def find_dofs(self, test=None, bc=None, boundary=True, dofrows=None, + check_vertices=True, check_facets=True, check_edges=True): + """Helper function for finding DOF indices for BC's. - Does not test for element interior DOFs since they are not typically included in - boundary conditions! Uses dofnum of 'u' variable. - - Remark: More convenient replacement for Assembler.dofnum_u.getdofs(). + Does not test for element interior DOFs since they are not typically + included in boundary conditions! Uses dofnum of 'u' variable. Parameters ---------- @@ -135,9 +133,6 @@ def essential_bc(self, test=None, bc=None, boundary=True, dofrows=None, check_ve I : np.array Set of DOF numbers set by the function """ - if self.mesh.dim() == 1: - raise Exception("Assembler.find_dofs not implemented for 1D mesh.") - if test is None: if self.mesh.dim() == 2: test = lambda x, y: 0*x + True @@ -194,15 +189,9 @@ def essential_bc(self, test=None, bc=None, boundary=True, dofrows=None, check_ve Fdofy = np.tile(my, (Fdofs.shape[0], 1)).flatten() locs = np.hstack((locs, np.vstack((Fdofx, Fdofy)))) else: - mx = 0.3333333*(self.mesh.p[0, self.mesh.facets[0, F]] + - self.mesh.p[0, self.mesh.facets[1, F]] + - self.mesh.p[0, self.mesh.facets[2, F]]) - my = 0.3333333*(self.mesh.p[1, self.mesh.facets[0, F]] + - self.mesh.p[1, self.mesh.facets[1, F]] + - self.mesh.p[1, self.mesh.facets[2, F]]) - mz = 0.3333333*(self.mesh.p[2, self.mesh.facets[0, F]] + - self.mesh.p[2, self.mesh.facets[1, F]] + - self.mesh.p[2, self.mesh.facets[2, F]]) + mx = np.sum(self.mesh.p[0, self.mesh.facets[:, F]], axis=0)/self.mesh.facets.shape[0] + my = np.sum(self.mesh.p[1, self.mesh.facets[:, F]], axis=0)/self.mesh.facets.shape[0] + mz = np.sum(self.mesh.p[2, self.mesh.facets[:, F]], axis=0)/self.mesh.facets.shape[0] Fdofx = np.tile(mx, (Fdofs.shape[0], 1)).flatten() Fdofy = np.tile(my, (Fdofs.shape[0], 1)).flatten() Fdofz = np.tile(mz, (Fdofs.shape[0], 1)).flatten() @@ -241,7 +230,8 @@ def essential_bc(self, test=None, bc=None, boundary=True, dofrows=None, check_ve elif self.mesh.dim() == 3: x[dofs] = bc(locs[0, :], locs[1, :], locs[2, :]) else: - raise NotImplementedError("Method essential_bc's not implemented for the given dimension.") + raise NotImplementedError("Method find_dofs not implemented " + + "for the given dimension.") return x, dofs diff --git a/skfem/mesh.py b/skfem/mesh.py index 441e2ff2a..26a963c66 100644 --- a/skfem/mesh.py +++ b/skfem/mesh.py @@ -13,9 +13,9 @@ >>> m.p.shape (2, 81) -Read a mesh generated using Gmsh. +Read a tetrahedral mesh generated using Gmsh. ->>> from skfem.extern_sfepy import read_gmsh +>>> from skfem.mesh_importers import read_gmsh >>> m = read_gmsh('examples/box.msh') >>> type(m) diff --git a/skfem/models/elasticity.py b/skfem/models/elasticity.py index 601363e64..30051688f 100644 --- a/skfem/models/elasticity.py +++ b/skfem/models/elasticity.py @@ -5,6 +5,30 @@ import numpy as np from skfem.utils import bilinear_form, linear_form +def basis_elasticity(m): + import skfem.mesh + from skfem.assembly import InteriorBasis + from skfem.element import ElementVectorH1, ElementHex1,\ + ElementTetP1, ElementTriP1 + from skfem.mapping import MappingAffine, MappingIsoparametric + if type(m) is skfem.mesh.MeshHex: + e1 = ElementHex1() + e = ElementVectorH1(e1) + map = MappingIsoparametric(m, e1) + return InteriorBasis(m, e, map, 3) + elif type(m) is skfem.mesh.MeshTet: + e1 = ElementTetP1() + e = ElementVectorH1(e1) + map = MappingAffine(m) + return InteriorBasis(m, e, map, 2) + elif type(m) is skfem.mesh.MeshTri: + e1 = ElementTriP1() + e = ElementVectorH1(e1) + map = MappingAffine(m) + return InteriorBasis(m, e, map, 2) + else: + raise NotImplementedError("basis_elasticity does not support the given mesh.") + def plane_strain(Lambda=1.0, Mu=1.0): @bilinear_form def weakform(u, du, v, dv, w): @@ -32,6 +56,47 @@ def Eps(dw): def lame_parameters(E, nu): return E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 - 2.0*nu)) +def boundary_prescribed(x=None, y=None, z=None): + """Supports only 1 DOF/node elements.""" + def prescribed(ind): + nonlocal x, y, z + def bc(basis): + nonlocal x, y, z + w = np.zeros(basis.dofnum.N) + D = np.array([]) + if x is not None: + if type(x) is int or type(x) is float: + X = lambda a, b, c: 0*a + x + else: + X = x + x1, D1 = basis.find_dofs(ind, bc=X, dofrows=[0], check_facets=False, check_edges=False) + w += x1 + D = np.concatenate((D, D1)) + if y is not None: + if type(y) == int or float: + Y = lambda a, b, c: 0*a + y + else: + Y = y + x2, D2 = basis.find_dofs(ind, bc=Y, dofrows=[1], check_facets=False, check_edges=False) + w += x2 + D = np.concatenate((D, D2)) + if z is not None: + if type(z) == int or float: + Z = lambda a, b, c: 0*a + z + else: + Z = z + x3, D3 = basis.find_dofs(ind, bc=Z, dofrows=[2], check_facets=False, check_edges=False) + w += x3 + D = np.concatenate((D, D3)) + return w, D + return bc + return prescribed + +def boundary_clamped(ind): + def bc(basis): + return basis.find_dofs(ind) + return bc + def linear_elasticity(Lambda=1.0, Mu=1.0): @bilinear_form def weakform(u, du, v, dv, w): diff --git a/skfem/utils.py b/skfem/utils.py index db7e7df6a..97aa986b3 100644 --- a/skfem/utils.py +++ b/skfem/utils.py @@ -85,7 +85,7 @@ def lin(v, dv, w): def smoothplot(x, basis, nref=3, m=None): - M, X = basis.refinterp(x, 3) + M, X = basis.refinterp(x, nref) if m is not None: ax = m.draw() M.plot(X, smooth=True, edgecolors='', ax=ax) @@ -210,3 +210,14 @@ def adaptive_theta(est, theta=0.5, max=None): else: return np.nonzero(theta*max < est)[0] +def initialize(basis, *bcs): + """Initialize a solution vector with boundary conditions in place.""" + y = np.zeros(basis.dofnum.N) + boundary_ix = np.array([]) + + for bc in bcs: + x, D = bc(basis) + y += x + boundary_ix = np.concatenate((boundary_ix, D)) + + return y, basis.dofnum.complement_dofs(boundary_ix)