3d elasticity example

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')

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())

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)
 

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

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))

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',

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
 

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)
 <class 'skfem.mesh.MeshTet'>

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):

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)