Multiscale finite element calculations in Python using SfePy
MMultiscale finite element calculations in Python usingSfePy
R. Cimrman a , V. Lukeˇs b , E. Rohan b a New Technologies - Research Centre, University of West Bohemia, Univerzitn´ı 8,301 00, Pilsen, Czech Republic b NTIS – New Technologies for the Information Society, Faculty of Applied Sciences,University of West Bohemia, Univerzitn´ı 8, 301 00, Pilsen, Czech Republic
Abstract
SfePy (Simple finite elements in Python) is a software for solving variouskinds of problems described by partial differential equations in one, two orthree spatial dimensions by the finite element method. Its source code ismostly (85%) Python and relies on fast vectorized operations provided bythe NumPy package. For a particular problem two interfaces can be used:a declarative application programming interface (API), where problem de-scription/definition files (Python modules) are used to define a calculation,and an imperative API, that can be used for interactive commands, or inscripts and libraries. After outlining the SfePy package development, thepaper introduces its implementation, structure and general features. Thecomponents for defining a partial differential equation are described using anexample of a simple heat conduction problem. Specifically, the declarativeAPI of SfePy is presented in the example. To illustrate one of SfePy’s mainassets, the framework for implementing complex multiscale models based onthe theory of homogenization, an example of a two-scale piezoelastic modelis presented, showing both the mathematical description of the problem andthe corresponding code.
Keywords: finite element method, multiscale simulations, piezoelasticity,SfePy, Python
Email addresses: [email protected] (R. Cimrman), [email protected] (V. Lukeˇs), [email protected] (E. Rohan)
Preprint submitted to Advances in Computational Mathematics October 2, 2018 a r X i v : . [ c s . M S ] O c t . Introduction SfePy ( http://sfepy.org ) is a software for solving systems of coupledpartial differential equations (PDEs) by the finite element method (FEM) in1D, 2D and 3D. It can be viewed both as a black-box PDE solver, and asa Python package which can be used for building custom applications. It isa multi-platform (Linux, Mac OS X, Windows) software released under theNew BSD license — the source code hosting, issue tracker and continuousintegration tools are available thanks to the GitHub development platform. SfePy has been employed by our group for a range of topics in mul-tiscale modelling in biomechanics and materials science, including multi-scale models of biological tissues (bone, muscle tissue with blood perfusion)[9, 38, 10, 34, 41, 40], a fish heart model with active contraction [27], com-putations of acoustic transmission coefficients across interfaces of arbitrarymicrostructure geometry [35], computations of phononic band gaps [39, 37],the finite element formulation of the Schroedinger equation [47, 12, 11], andother applications.In Section 2 the programming language choice is discussed and the
SfePy project development is described. In Section 3 an overview of the pack-age is given, and a simple example definition (a heat conduction problem)is presented. Section 4 introduces the
SfePy ’s homogenization engine —a sub-package for defining complex multiscale problems using a simple do-main specific language, based on Python’s dictionaries. Finally, in Section 5a complex example — a multiscale numerical simulation of a piezoelectricstructure — is shown as expressed in the homogenization engine syntax.
2. Development
The code has been written primarily in the Python programming lan-guage . Python is a very high-level interpreted programming language, thathas a number of features appealing to scientists (non-IT), such as: a clean,easy-to-read syntax, no manual memory management, a huge standard li-brary, a very good interoperability with other languages (C, fortran), anda large and friendly scientific computing community. It allows both fastexploration of various ideas and efficient implementation, thanks to many SfePy sources (version 2018.3) GitHub statistics: Python 85.1%, C 14.6%, other: 0.3%.
SfePy project uses Git [20] for source code management and GitHubdevelopment platform [21] for the source code hosting and developer interac-tion, see https://github.com/sfepy/sfepy , similarly to many other scien-tific Python tools. Travis CI [46] is used (via GitHub) for running automatictests after every uploaded commit. The developers and users of the softwarecan communicate using the mailing list “[email protected]”. The source codeand the package usage and development are documented with the help of theSphinx documentation generator [45] that is also used to generate the pagesof the
SfePy project web-site.The version 2018.3 has been released in September 17, 2018, and itsgit-controlled sources contained 897 total files, 721447 total lines (1539277added, 817830 removed) and 6386 commits done by 24 authors . About120000 lines (16 %) are the source code, the other lines belong to the finiteelement meshes, documentation etc.
3. Description
In this section we briefly outline the package implementation, structureand general features. The components for defining a PDE to solve are de-scribed using a simple example in Section 3.4.
Because Python is an interpreted language, and the standard implemen-tation (CPython) has slow loops, several approaches are used in the codeto achieve good (C-like) performance. For speed in general, it relies on fastvectorized operations provided by NumPy arrays [30], with significant use ofadvanced features such as broadcasting. C and Cython [6] are used in placeswhere vectorization is not possible, or is too difficult/unreadable. Most of the authors contributed only one or a few commits. fePy relies on a number of packages of the scientific Python softwarestack, namely: SciPy [25] for sparse matrices, solvers and algorithms, Mat-plotlib [24] for 2D plots, Mayavi [32] for 3D plots and simple post-processingGUI, PyTables [31] for HDF5 file format support, SymPy [29] for symbolicoperations/code generation, igakit (a part of PetIGA [16]) for working withNURBS bases of the isogeometric analysis [14, 5] etc.Besides the vectorized operations of NumPy, other performance gainsare enabled by using very efficient solvers with a Python interface, such asUMFPACK [18] + scikit-umfpack [42], MUMPS [3] or PETSc [4]. SfePy canrun in parallel, using PETSc + petsc4py [17] and mpi4py [15] packages. Weemploy the separation of concerns strategy with respect to parallelism: mostof the
SfePy code is serial, and there is a single dedicated module in
SfePy ,that, together with the flexible way of computing weak form integrals on anysubdomain (see below), allows parallel assembling of the discrete systemsand their subsequent solution. In SfePy , the equations are not given in a fully symbolic way, as in, forexample, Fenics or Firedrake projects , but a simpler approach is used: the SfePy package comes with a database of predefined terms. A term is thesmallest unit that can be used to build equations . It corresponds to a weakformulation integral over a (sub)domain and takes usually several arguments:(optional) material parameters, a single virtual (or test) function variable andzero or more state (or unknown) variables. The available terms are listed atour web site [43], currently there are 118 terms.The high-level code that handles a PDE discretization in
SfePy is inde-pendent of a method of domain or variable discretization. For each particularmethod of discretization, there is a sub-package that implements the spe-cific functionality (degrees of freedom management, selection of subdomains,reference domain mappings, etc.). This abstraction allows adding variousdiscretization methods. The following ones are currently implemented: • the finite element method on 1D line, 2D area (triangle, rectangle) and3D volume (tetrahedron, hexahedron) finite elements; with two kindsof polynomial bases: Both use the Unified Form Language from Fenics. the classical nodal (Lagrange) basis that can be used with allsupported element/cell types; – the hierarchical (Lobatto) basis [44] that can be used with tensor-product elements (line, rectangle, hexahedron).The basis function polynomials of an arbitrary order (theoretically, seelimitations below) as well as the corresponding quadrature rules aresupported. The Lagrange basis is implemented in C/Cython, while theLobatto basis functions use a C code generated using SymPy. • the isogeometric analysis [14] with a NURBS or B-spline basis, imple-mented using the B´ezier extraction approach [5], currently limited tosingle NURBS patch domains, see [7].All the basis functions listed above support the H function spaces only( (cid:126)H (curl), (cid:126)H (div) spaces are not currently implemented). In addition tothe above elements, two structural elements are implemented (using SfePy terms): the hyperelastic Mooney-Rivlin membrane [48], and the shell10xelement term based on the Reissner-Mindlin theory [49]. can solve many problems described by PDEs in the weak form. Fora particular problem, there are two interfaces that can be used: • a declarative API, where problem description/definition files (Pythonmodules) are used to define a calculation; • an imperative API, that can be used for interactive commands, or inscripts and libraries.Both the above APIs closely correspond to the mathematical description ofthe weak form PDEs. An advanced use of the declarative API is demon-strated in the piezoelectric model example in Section 5.The declarative API involves almost no programming besides using basicPython data types (dicts, lists, tuples, strings, etc.) and allows a lazy defini-tion of the problem, called problem configuration , as well as a manipulationwith the problem configuration. Prior to a problem solution, the problemconfiguration is automatically translated into a problem object using theimperative API. The SfePy package contains several top-level scripts thatcan be used to run simulations defined using the declarative API. The twocommon ones are: 5 the simple.py script that allows running regular calculations of PDEs, • the homogen.py script that allows running the homogenization engineto compute effective material parameters, see Section 4.The imperative API allows immediate evaluation of expressions, and thussupports interactive exploration or inspection of the FE data. It is also morepowerful than the declarative API as a user is free to perform non-predefinedtasks. The problems defined using the imperative API usually have a main() function and can be run directly using the Python interpreter.In the both cases, a problem definition is a Python module, so all thepower of Python (and supporting SfePy modules) is available when neededfor complex problems.
Systems of PDEs are defined using keywords or classes corresponding tomathematical objects present in the weak formulation of the PDEs. Here weillustrate the components of the problem definition using a simple example.We wish to solve a heat conduction problem, that can be written in the weakform as follows: Find the temperature u ∈ H (Ω) such that for all v ∈ H (Ω)holds (cid:90) Ω v ∂u∂t + (cid:90) Ω c ∇ v · ∇ u = 0 , ∀ v , u ( x,
0) = g ( x ) , u ( x, t ) = (cid:26) − x ∈ Γ left , x ∈ Γ right , where c is a material parameter (a thermal diffusivity). Below we show thedeclarative way of defining the ingredients necessary to solve this problemwith SfePy . For complete examples illustrating both the declarative andimperative APIs, see the accompanying dataset [8]. • The domain Ω has to be discretized, resulting in a finite element mesh .The mesh can be loaded from a file (generated by external tools) orgenerated by the code (simple shapes). filename_mesh = ’meshes/3d/cylinder.mesh’ • Regions serve as domains of integration and allow defining boundaryand initial conditions. Subdomains of various topological dimensioncan be defined. The mesh/domain and region handling uses a C datastructure adapted from [28]. The following code defines the domain Ωand the boundaries Γ left , Γ right . 6 egions = {’Omega’ : ’all’,’Left’ : (’vertices in (x < 0.00001)’, ’facet’),’Right’ : (’vertices in (x > 0.099999)’, ’facet’),} • Fields correspond to the discrete function spaces and are defined usingthe (numerical data type, number of components, region name, approx-imation order) tuple. A field can be defined on the whole domain, ona volume (cell) subdomain or on a surface (facet) region. fields = {’temperature’ : (’real’, 1, ’Omega’, 1),} • The fields (FE spaces) can be used to define variables . Variablescome in three flavors: unknown field for state variables, test field for test (virtual) variables and parameter field for variables withknown values of degrees of freedom (DOFs). The definition items foran unknown variable definition are: ( ’unknown field’ , field name,order in global vector, [optional history size]) . In the snippet below,the history size is 1, as the previous time step state is required for thenumerical time derivative. For a test variable, the last item is the nameof the corresponding unknown variable. variables = {’u’ : (’unknown field’, ’temperature’, 0, 1),’v’ : (’test field’, ’temperature’, ’u’),} • Materials correspond to all parameters defined point-wise in quadra-ture points, that can be given either as constants, or as general func-tions of time and quadrature point coordinates. Here we just definethe constant parameter c , as a part of the material ’m’ . materials = {’m’ : ({’c’ : 1.0e-5},),} Similarly to materials, the Dirichlet (essential) boundary conditions can be defined using constants or general functions of time and coordi-nates. In our case we set the values of u to 2 and -2 on Γ left and Γ right ,respectively. ebcs = {’u1’ : (’Left’, {’u.0’ : 2.0}),’u2’ : (’Right’, {’u.0’ : -2.0}),} • The initial conditions can be defined analogously, here we illustratehow to use a function. The conditions are applied in the whole domainΩ. The code assumes NumPy was imported ( import numpy as np ),and ic max is a constant defined outside the function. def get_ic(coors, ic):x, y, z = coors.Treturn 2 - 40.0 * x + ic_max * np.sin(4 * np.pi * x / 0.1)functions = {’get_ic’ : (get_ic,),}ics = {’ic’ : (’Omega’, {’u.0’ : ’get_ic’}),} • The PDEs can be built as a linear combination of many predefinedterms. Each term has its quadrature order and its region of integration.The integral specifies a numerical quadrature order. integrals = {’i’ : 2,}equations = {’Temperature’ : """dw_volume_dot.i.Omega(v, du/dt)+ dw_laplace.i.Omega(m.c, v, u) = 0"""} • A complete example contains additional information (notably a config-uration of solvers), see [8]. There, the above code snippets are used in8he heat cond declarative.py file. The simulation is then launchedusing the python
Figure 1: Initial and final snapshots of the temperature evolution. has no meshing capabilities besides several simple mesh genera-tors (a block mesh generator, an open/closed cylinder mesh generator, amesh generator from CT image data), but several operations like mergingof matching meshes are supported. The FE mesh needs to be provided ina file in one of the supported formats, notably the legacy VTK format [26].The results are stored in legacy VTK files, or in, usually in case of time-dependent problems with many time steps, custom HDF5 [22] files. Manystandard open-source tools can be used to display the VTK files, namelyParaview [23], or Mayavi [32], see data workflow of a simulation in Fig. 2.Mayavi is supported directly within
SfePy via the postproc.py script. The extractor.py script is provided to extract/convert the HDF5 file content toVTK. provides and uses a unified interface to many standard codes, forexample UMFPACK [18], MUMPS [3], PETSc [4], Pysparse [19] as well asthe solvers available in SciPy. Various solver classes are supported: time-stepping, nonlinear, linear, eigenvalue problem and optimization solvers. Anautomatically generated list of all the supported solvers can be found at the
SfePy web site [43]. 9
E solver problem descriptionfile (.py)FE mesh file (.vtk, .msh, .hmascii,.h5, .bdf, .neu, ...) MayaviParaview...NetgenHypermesh ...
GmshANSYS post-processingpre-processing output file(.vtk or .h5)
Figure 2: Data exchange between
SfePy and external pre-processing and post-processingtools.
Besides external solvers, several solvers are implemented directly in
SfePy ,for example: • ts.simple : Implicit time stepping solver with a fixed time step, suit-able also for quasistatic problems. • ts.newmark , ts.bathe , ts.generalized alpha : Solve elastodynamicsproblems by the Newmark, Bathe, generalized- α methods, respectively. • nls.newton : The Newton nonlinear solver with a backtracking line-search.A typical problem solution, when using the declarative problem definitionAPI, then involves calling a time-stepping solver that calls a nonlinear solverin each time step, which, in turn, calls a linear solver in the nonlinear iter-ations. A unified approach is used here: for stationary problems, a dummytime-stepping solver ( ts.stationary ) is automatically created. Similarly,a nonlinear solver is used to solve both the linear and non-linear problems.This simplifies imposing non-homogeneous Dirichlet boundary conditions: inthe context of a nonlinear solver, the increment of the constrained DOFs isalways zero and the non-zero boundary values can be ignored during theassembling. 10 .7. Limitations The limitations can be split into two groups. The first group is relatedto limited number of developers and our research focus: certain features aremissing, because they do not fall into our field of research (e.g. the vectorfinite elements). The limitations in the second group are more fundamental.Because the code relies on vectorization provided by NumPy, the code triesto work on all cells in a region in each operation: for example, all localfinite element matrices are evaluated in a vectorized way into a single largeNumPy array, and then assembled to a SciPy’s sparse matrix. This placesa restriction on a practically usable order of the basis function polynomials,especially for 3D hexahedrons, where orders greater than 4 are not practicallyusable. Using NumPy’s arrays places another restriction: data homogeneity.So, for examples, the FE basis polynomial order has to be uniform overthe whole (sub)domain where a field is defined. This is incompatible withan adaptive mesh refinement, especially the hp -adaptivity. Note that SfePy supports meshes with level-1 hanging nodes, so a limited h -adaptivity ispartly possible.
4. Homogenization engine
In this section we briefly outline the approach to solving multiscale prob-lems based on the theory of homogenization [1, 13]. The main asset of thehomogenization approach is that a homogenized model can take into ac-count various details at the microstructure scale (topology, heterogeneousmaterial parameters, etc.) without actually meshing those detailed featureson a macroscopic domain, which would lead to an extremely large problem.The homogenization engine is a feature that is in our opinion unique to
SfePy . It has been developed to allow an easy and flexible formulation ofproblems arising from use of the homogenization theory applied to stronglyheterogeneous multiscale material models. Such models, as can be seen inthe non-trivial example in Section 5, can have a complicated data flow anddependencies of various subproblems involved in the definition.A typical homogenization procedure for a linear problem involves thefollowing steps: The situation is much more complicated for nonlinear problems: a microproblem needsto be typically solved in every macroscopic integration point. This mode is also supportedin
SfePy .
11. Compute characteristic (corrector) functions by solving auxiliary cor-rector problems on a reference periodic cell domain that describes themicrostructure (exactly or in a statistical sense).2. Using the corrector functions, evaluate the homogenized coefficients.Those coefficients correspond to effective macroscopic properties of thematerial with the given microstructure as the microstructure charac-teristic scale tends to zero.3. Solve the homogenized model with the obtained effective homogenizedcoefficients on a macroscopic domain.There can be a number of the auxiliary corrector problems as well as thehomogenized coefficients. The homogenization engine allows to describe theirrelationships and the dependencies among them and resolves the problems ina correct order automatically. A complete problem is described in one or moreproblem definition files using the declarative API, and data dependencies aredescribed using Python dictionaries as a small domain-specific language.The homogenization engine allows to solve microscopic subproblems andevaluate homogenized coefficients in parallel with use of either multithread-ing or multiprocessing features of a computer system. The distribution ofmicroproblems between multiple threads or CPUs is governed by a functionthat puts the microproblems with resolved dependencies into the work queue,collects solved microscopic solutions and updates a dependency table accord-ing to the obtained results. Workers, i.e. threads or CPU cores, solve tasksfrom the work queue until it is empty. If the same microproblem needs to besolved multiple times with different parameters, typically for nonlinear prob-lems, the total amount of microscopic tasks is divided into several chunksthat are distributed to multiple workers. In the case of multiprocessing, theMPI library is used to communicate between computational nodes.The parallel computation is crucial for nonlinear problems where micro-problems have to be resolved in all macroscopic integration points in manytime or iteration steps.
5. Multiscale numerical simulation of piezoelectric structure
The complex multiscale model, solved by the means of the homogeniza-tion method, and its implementation in
SfePy are presented in this section.12 .1. Mathematical model of piezoelectric media
We consider a porous piezoelectric medium which consists of a piezoelec-tric matrix, embedded metallic electrodes (conductors) and void inclusions.These components are arranged in a periodic lattice so that the medium canbe generated by copies of the reference unit cell, see Fig.3. The mechanicalbehavior of such a structure can be described using the two-scale asymptotichomogenization method, see [1, 13]. The quantities oscillating within theheterogeneous structure with the period equal to the size of the periodic unitare labelled by the superscript ε in the subsequent text.The mechanical properties of the piezoelectric solid are given by the fol-lowing constitutive equations σ εij ( (cid:126)u ε , ϕ ε ) = A εijkl e kl ( (cid:126)u ε ) − g εkij E k ( ϕ ε ) ,D εk ( (cid:126)u ε , ϕ ε ) = g εkij e ij ( (cid:126)u ε ) + d εkl E l ( ϕ ε ) , (1)which express the dependencies of the Cauchy stress tensor σ ε and the electricdisplacement (cid:126)D ε on the strain tensor (cid:126)e ( (cid:126)u ε ) = (cid:0) ∇ (cid:126)u ε + ( ∇ (cid:126)u ε ) T (cid:1) , where (cid:126)u ε isthe displacement field, and on the electric field (cid:126)E ( ϕ ε ) = ∇ ϕ ε , where ϕ ε isthe electric potential. On the right hand side of (1), we have the fourth-orderelastic tensor A εijkl ( A εijkl = A εklij = A εjilk ), the third-order tensor g εkij ( g εkij = g εkji ), which couples mechanical and electric quantities, and the permeabilitytensor d εkl .The quasi-static problem of the piezoelectric medium is given by thefollowing equilibrium equations −∇ · σ ε ( (cid:126)u ε , ϕ ε ) = (cid:126)f ε , in Ω εmc , −∇ · (cid:126)D ε ( (cid:126)u ε , ϕ ε ) = q εE , in Ω εm , (2)and by the boundary conditions (cid:126)n · σ ε = (cid:126)h ε on Γ εσ , (cid:126)n · (cid:126)D ε = (cid:37) εE on Γ ε(cid:126)D ,(cid:126)u ε = ¯ (cid:126)u on Γ ε(cid:126)u , ϕ ε = ¯ ϕ on Γ εϕ , (3)where (cid:126)f ε , (cid:126)h ε are the volume and surface forces, q εE , (cid:37) εE are the volume andsurface charges and ¯ (cid:126)u , ¯ ϕ are the prescribed displacements and electric po-tential, respectively. The piezo-elastic medium occupies an open boundedregion Ω ∈ R which is decomposed into several non-overlapping parts: thepiezoelectric elastic matrix Ω εm , conductive elastic parts Ω εc = (cid:83) k Ω k,εc and13 igure 3: The scheme of the representative periodic cell decomposition and the generatedperiodic structure. isolated void inclusions Ω εo , see Fig. 3. The elastic part, i.e. the matrix andthe conductors, is denoted by Ω εmc = Ω εm ∪ Ω εc .The weak formulation of the problem stated above can be written as:Given volume and surface forces (cid:126)f ε , (cid:126)h ε and volume and surface charges q εE , (cid:37) εE , find (cid:126)u ε ∈ U (Ω εmc ), ϕ ε ∈ V (Ω εm ) such that for all (cid:126)v ∈ U (Ω εmc ), ψ ∈ V (Ω εm ) (cid:90) Ω εm [ (cid:126)A ε (cid:126)e ( (cid:126)u ε ) − ( (cid:126)g ε ) T · ∇ ϕ ε ] : (cid:126)e ( (cid:126)v ) d V + (cid:90) Ω εc [ (cid:126)A ε (cid:126)e ( (cid:126)u ε )] : (cid:126)e ( (cid:126)v ) d V = (cid:90) Γ εσ (cid:126)h · (cid:126)v d S + (cid:90) Ω εmc (cid:126)f ε · (cid:126)v d V , (cid:90) Ω εm [ (cid:126)g ε : (cid:126)e ( (cid:126)u ε ) + (cid:126)d ε · ∇ ϕ ε ] · ∇ ψ d V = (cid:90) Γ ε(cid:126)D (cid:37) εE ψ d S + (cid:90) Ω εm q εE ψ d V . (4)Symbols U , U , V , V denote admissibility sets, where U , V are the setswith zero trace on the Dirichlet boundary. Further details can be found in[36]. We apply the standard homogenization techniques, cf. [1] or [13], to theproblem (4). It results in the limit model for ε −→
0, where ε is the scaleparameter relating the microscopic and macroscopic length scales. The ho-mogenization process leads to local microscopic problems, defined within areference periodic cell, and to the global problem describing the behaviorof the homogenized medium at the macroscopic level. The global probleminvolves the homogenized material coefficients which are evaluated using thesolutions of the local problems. Due to linearity of the problem, the micro-scopic and macroscopic problems are decoupled.14s we assume given potentials ¯ ϕ k in each of the electrode networks, thedielectric properties must be appropriately rescaled in order to preserve thefinite electric field for the limit ε −→ (cid:126)g ε = ε ¯ (cid:126)g , (cid:126)d ε = ε ¯ (cid:126)d , cf. [36]. Local problems and homogenized coeffficients.
The local microscopic responsesof the piezoelectric structure are given by the following sub-problems whichare solved within the periodic reference cell Y , see Fig. 3, that is decomposedsimilarly to the decomposition of domain Ω: • Find ω ij ∈ (cid:126)H ( Y mc ), η ij ∈ H ( Y m ) such that for all (cid:126)v ∈ (cid:126)H ( Y mc ), ψ ∈ H ( Y m ) and for any i, j = 1 , , (cid:90) Y mc (cid:104) (cid:126)A(cid:126)e (cid:0) ω ij + Π ij (cid:1)(cid:105) : (cid:126)e ( (cid:126)v ) d V − (cid:90) Y m (cid:2) ¯ (cid:126)g T · ∇ η ij (cid:3) : (cid:126)e ( (cid:126)v ) d V = 0 , (cid:90) Y m (cid:104) ¯ (cid:126)g : (cid:126)e (cid:0) ω ij + Π ij (cid:1) + ¯ (cid:126)d · ∇ η ij (cid:105) · ∇ ψ d V = 0 , (5)where Π ijk = y j δ ik . • Find ˆ ω k ∈ (cid:126)H ( Y mc ), ˆ η k ∈ H ,k ( Y m ) such that for all (cid:126)v ∈ (cid:126)H ( Y mc ), ψ ∈ H ( Y m ) and for any k = 1 , , . . . , k c ( k c is the number of conductors) (cid:90) Y mc (cid:104) (cid:126)A(cid:126)e (cid:0) ˆ ω k (cid:1)(cid:105) : (cid:126)e ( (cid:126)v ) d V − (cid:90) Y m (cid:2) ¯ (cid:126)g T · ∇ ˆ η k (cid:3) : (cid:126)e ( (cid:126)v ) d V = 0 , (cid:90) Y m (cid:104) ¯ (cid:126)g : (cid:126)e (cid:0) ˆ ω k (cid:1) + ¯ (cid:126)d · ∇ ˆ η k (cid:105) · ∇ ψ d V = 0 . (6)The microscopic sub-problems are solved with the periodic boundary condi-tions and ˆ η k = δ ki on Γ imc for i = 1 , , . . . , k c , Γ imc = Y m ∩ Y ic is the interfacebetween the matrix part Y m and i -th conductor Y ic . By (cid:126)H we refer to theSobolev space of Y -periodic functions, H ,k reflects the above mentionedinterface condition on Γ imc and H is the set of functions which are equal tozero on Γ mc .With the characteristic responses ω ij , η ij and ˆ ω k , ˆ η k obtained by solv-ing (5) and (6), the homogenized material coefficients (cid:126)A H and (cid:126)P H,k can be15valuated using the following expressions: A Hijkl = 1 | Y | (cid:20)(cid:90) Y mc (cid:104) (cid:126)A(cid:126)e (cid:0) ω ij + Π ij (cid:1)(cid:105) : (cid:126)e (cid:0) ω kl + Π kl (cid:1) d V + (cid:90) Y m ¯ (cid:126)d ∇ η ij · ∇ η kl d V (cid:21) ,P H,kij = 1 | Y | (cid:20)(cid:90) Y mc (cid:104) (cid:126)A(cid:126)e (cid:0) ˆ ω k (cid:1)(cid:105) : (cid:126)e (cid:0) Π ij (cid:1) d V − (cid:90) Y m (cid:2) ¯ (cid:126)g : (cid:126)e (cid:0) Π ij (cid:1)(cid:3) · ∇ ˆ η k d V (cid:21) . (7) Macroscopic problem.
The global macroscopic problem is defined in terms ofthe homogenized coefficients as: Find the macroscopic displacements (cid:126)u ∈ U (Ω) such that for all (cid:126)v ∈ U (Ω) (cid:90) Ω [ (cid:126)A H (cid:126)e (cid:0) (cid:126)u (cid:1) ] : (cid:126)e ( (cid:126)v ) d V = − (cid:90) Ω (cid:126)e ( (cid:126)v ) : (cid:88) k (cid:126)P H,k ¯ ϕ k d V . (8)We assumed that (cid:37) E = 0, otherwise we would need an extra coefficient relatedto the surface charge, see [36]. In this section we illustrate the use of
SfePy ’s homogenization engine inthe following setting. The macroscopic problem described by (8) is solved inthe domain Ω, depicted in Fig. 4 right, that is fixed at its left face ( u i = 0 for i = 1 , , left ). No volume and surface forces or charges are applied andthe deformation of the macroscopic sample is invoked only by the prescribedelectrical potential ¯ ϕ = ± V in the two embedded conductor networks.The geometry of the representative volume element, which is used to solvethe microscopic problems (5), (6), is depicted in Fig. 4 left. The materialparameters of the piezoelectric elastic matrix, made of barium-titanite, andmetallic conductors are summarized in Table 1.The results of the multiscale numerical simulation are shown in Figs. 5, 6.The macroscopic strain field and the deformed macroscopic sample (defor-mation scaled by factor 100) are presented in Fig. 5. Using the macroscopicsolution and the characteristic responses, we can reconstruct the fields at themicroscopic level for a given finite size ε in a chosen part of the macroscopicdomain. The reconstructed strain field and the deformed microstructure(deformation scaled by factor 10) are shown in Fig. 6 left, the reconstructedelectric field is depicted in Fig. 6 right. See [36] for comparison of accuracy ofthe reconstructed solutions with a full, much more demanding, simulation. In16 iezoelectric matrix Elasticity – transverse isotropy [GPa]: A A A A A A ( A = A , A = A , 15.040 14.550 6.560 6.590 4.240 4.390 A = A )Piezo-coupling [C/m ]: g g g g ( g = g , g = g ) -4.322 -4.322 17.360 11.404Dielectricity [10 − C/Vm]: d d ( d = d ) 1.284 1.505 Metallic conductors
Elasticity – linear isotropy: E [GPa] ν [-]200.0 0.25 Table 1: Properties of the piezoelectric matrix and metallic conductors.Figure 4: Left: the geometry of the reference periodic cell Y ; Right: the macroscopicdomain Ω. the above article, the full (reference) simulation has approximately 4 . × degrees of freedom while the microscopic problem has only 741 DOFs andthe macroscopic problem 577 DOFs. The linear multiscale analysis defined above is performed in
SfePy in twosteps:1. The local microscopic sub-problems (5), (6) are solved using the ho-mogenization engine, see Section 4. The engine is also used to evaluatethe homogenized coefficients according to (7).2. The global macroscopic problem (8) is solved, the known homogenizedcoefficients are employed. 17 igure 5: The deformed macroscopic sample (deformation scaled by factor 100) and themagnitude of macroscopic strain field.Figure 6: Left: the deformed microscopic structure (deformation scaled by factor 10) andthe magnitude of reconstructed strain field; Right: the magnitude of the reconstructedelectric field. equations = {’balance_of_forces’: """dw_lin_elastic.i2.Omega(hom.A, v, u)= - dw_lin_prestress.i2.Omega(hom.Pf, v)""",} where hom.A stands for the homogenized coefficients (cid:126)A H and hom.Pf is equalto (cid:80) k (cid:126)P H,k ¯ ϕ k , i.e. the sum of the coefficients (cid:126)P H,k multiplied by the pre-scribed electrical potentials ¯ ϕ k . The homogenized material hom is declaredas a function, which calls the homogenization engine (via SfePy ’s built-infunction get homog coefs linear() in the following code) and returns thecalculated homogenized parameters. In the case of a linear problem, the samevalues are valid in all quadrature points of a macroscopic domain ( coors ar-gument in the function below). def get_homog_fun(ts, coors, mode, **kwargs):...coefs = get_homog_coefs_linear(0, 0, None,micro_filename=’piezo_elasticity_micro.py’,coefs_filename=’coefs_piezo.h5’)Pf = coefs[’P1’] * bar_phi[0] + coefs[’P2’] * bar_phi[1]nqp = coors.shape[0]out = {’A’: np.tile(coefs[’A’], (nqp, 1, 1)),’Pf’: np.tile(Pf[:, np.newaxis], (nqp, 1, 1)),}return outfunctions = {’get_homog’: (get_homog_fun,),}materials = {’hom’: ’get_homog’,}
To define the microscopic sub-problems which are solved by the homog-enization engine the following fields and variables are needed:19 ields = {’displacement’: (’real’, ’vector’, ’Ymc’, 1),’potential’: (’real’, ’scalar’, ’Ym’, 1),}variables = { ’u’: (’unknown field’, ’displacement’),’v’: (’test field’, ’displacement’, ’u’),’Pi_u’: (’parameter field’, ’displacement’, ’u’),’U1’: (’parameter field’, ’displacement’, ’(set-to-None)’),’U2’: (’parameter field’, ’displacement’, ’(set-to-None)’), ’r’: (’unknown field’, ’potential’),’s’: (’test field’, ’potential’, ’r’),’Pi_r’: (’parameter field’, ’potential’, ’r’),’R1’: (’parameter field’, ’potential’, ’(set-to-None)’),’R2’: (’parameter field’, ’potential’, ’(set-to-None)’),}
The material parameters of the elastic matrix ( Y m ) and the metalic conduc-tors ( Y c ) are defined as follows, see Table 1: materials = {’elastic’: ({’D’: {’Ym’: np.array([[1.504, 0.656, 0.659, 0, 0, 0],[0.656, 1.504, 0.659, 0, 0, 0],[0.659, 0.659, 1.455, 0, 0, 0],[0, 0, 0, 0.424, 0, 0],[0, 0, 0, 0, 0.439, 0],[0, 0, 0, 0, 0, 0.439]]) * 1e11,’Yc’: stiffness_from_youngpoisson(3, 200e9, 0.25)}},),’piezo’: ({’g’: np.array([[0, 0, 0, 0, 11.404, 0],[0, 0, 0, 0, 0, 11.404],[-4.322, -4.322, 17.360, 0, 0, 0]]),’d’: np.array([[1.284, 0, 0],[0, 1.284, 0],
0, 0, 1.505]]) * 1e-8},),}
The homogenized coefficients (cid:126)A H can be introduced as import sfepy.homogenization.coefs_base as cbcoefs = {’A1’: {’requires’: [’pis_u’, ’omega_ij’],’expression’: ’dw_lin_elastic.i2.Ymc(elastic.D, U1, U2)’,’set_variables’: [(’U1’, (’omega_ij’, ’pis_u’), ’u’),(’U2’, (’omega_ij’, ’pis_u’), ’u’)],’class’: cb.CoefSymSym,},’A2’: {’requires’: [’omega_ij’],’expression’: ’dw_diffusion.i2.Ym(piezo.d, R1, R2)’,’set_variables’: [(’R1’, ’omega_ij’, ’r’),(’R2’, ’omega_ij’, ’r’)],’class’: cb.CoefSymSym,},’A’: {’requires’: [’c.A1’, ’c.A2’],’expression’: ’c.A1 + c.A2’,’class’: cb.CoefEval,},} where we follow the expression (7 ) which consists of two integrals over do-mains Y mc (matrix + conductors) and Y c (conductors). The definition ofeach coefficient has these parts: requires – the names of correctors neededfor evaluation, expression – the expression to be evaluated, class – thecoefficient class; it determines the way of evaluation and the resulting ma-trix/array shape. In our case, the class of A1 and A2 is CoefSymSym : it meansthat the resulting coefficients are the fourth-order tensors in the symetricstorage, e.g. sym × sym matrices, where sym is the number of componentsin a symmetric stress/strain vector. Class CoefEval is used to evaluatea simple mathematical expression, in our example, the summation of A1 , A2 . In set variables section we say how to substitute the correctors into21he variables employed in the expression. For example, the code (’U1’,(’omega ij’, ’pis u’), ’u’) is interpretted as: U ω K + Π K , where ω K is stored in omega ij[’u’] , Π K in pis u[’u’] and K is the multi-index attaining 11 , , , , ,
23 for a 3D problem because of the used
CoefSymSym class. In a similar way, the coefficients (cid:126)P H, can be introducedas coefs.update({’P1_1’: {’requires’: [’pis_u’, ’omega_k1’],’expression’: ’dw_lin_elastic.i2.Ymc(elastic.D, U1, U2)’,’set_variables’: [(’U1’, ’omega_k1’, ’u’),(’U2’, ’pis_u’, ’u’)],’class’: cb.CoefSym,},’P1_2’: {’requires’: [’pis_u’, ’omega_k1’],’expression’: ’dw_piezo_coupling.i2.Ym(piezo.g, U1, R1)’,’set_variables’: [(’R1’, ’omega_k1’, ’r’),(’U1’, ’pis_u’, ’u’)],’class’: cb.CoefSym,},’P1’: {’requires’: [’c.P1_1’, ’c.P1_2’],’expression’: ’c.P1_1 - c.P1_2’,’class’: cb.CoefEval,},}) Here, the
CoefSym class is employed due to the second-order coefficient whichcan be represented as a vector with dimension sym .The required correctors omega ij , see (5), and pis u are defined as fol-lows: requirements = {’pis_u’: {’variables’: [’u’],’class’: cb.ShapeDimDim,}, omega_ij’: {’requires’: [’pis_u’],’ebcs’: [’fixed_u’, ’fixed_r’],’epbcs’: [’p_ux’, ’p_uy’, ’p_uz’, ’p_rx’, ’p_ry’, ’p_rz’],’equations’: {’eq1’: """dw_lin_elastic.i2.Ymc(elastic.D, v, u)- dw_piezo_coupling.i2.Ym(piezo.g, v, r)= -dw_lin_elastic.i2.Ymc(elastic.D, v, Pi_u)""",’eq2’: """dw_piezo_coupling.i2.Ym(piezo.g, u, s)+ dw_diffusion.i2.Ym(piezo.d, s, r)= -dw_piezo_coupling.i2.Ym(piezo.g, Pi_u, s)""",},’set_variables’: [(’Pi_u’, ’pis_u’, ’u’)],’class’: cb.CorrDimDim,’dump_variables’: [’u’, ’r’],},} The class
ShapeDimDim is used to define the symbol Π ijk = y j δ ik and CorrDimDim ensures the corrector with dim × dim components, dim is the space dimen-sion. Note that a corrector can also depend on another corrector as in thecode above, where pis u is required to solve omega ij . The correctors intro-duced in (6) can be defined in SfePy as requirements.update({’omega_k1: {’requires’: [’pis_r’],’ebcs’: [’fixed_u’, ’fixed_interface’],’epbcs’: [’p_ux’, ’p_uy’, ’p_uz’, ’p_rx’, ’p_ry’, ’p_rz’],’equations’: {’eq1’: """dw_lin_elastic.i2.Ymc(elastic.D, v, u)- dw_piezo_coupling.i2.Ym(piezo.g, v, r) = 0""",’eq2’: """dw_piezo_coupling.i2.Ym(piezo.g, u, s)+ dw_diffusion.i2.Ym(piezo.d, s, r) = 0"""},’class’: cb.CorrOne,’dump_variables’: [’u’, ’r’],},}) CoefOne corresponds to the scalar corrector function. The correc-tors are solved with the periodic boundary conditions defined in the lineswith the keyword epbcs and the Dirichlet (essential) boundary conditionsdefined in the lines with the keyword ebcs .The multiscale simulation can be run by calling the simple.py script,see Section 3.5, with the name of the description file for the macroscopicproblem as a script parameter. The script runs the simulation at the macro-scopic level and invokes the homogenization engine through the materialfunction. The full sources of this example can be found in the
SfePy pack-age in examples/multiphysics/ : piezo elasticity macro.py defines themacroscopic problem, and piezo elasticity micro.py defines the compu-tations on the reference periodic cell of the microstructure. For the versionof the sources used in this article see [8].
6. Conclusion
We introduced the open source finite element package
SfePy , a code writ-ten (mostly) in Python for solving various kinds of problems described bypartial differential equations and discretized by the finite element method.The design of the code was discussed and illustrated using a simple heatconduction example.Special attention was devoted to the description of the
SfePy ’s homog-enization engine, a sub-package for defining complex multiscale problems.This feature was introduced in a tutorial-like form using a multiscale numer-ical simulation of a piezoelectric structure.For the complete code of the examples presented, together with the re-quired FE meshes and the 2018.3 version of
SfePy , see [8]. Further documen-tation and many more examples of
SfePy use can be found on the project’sweb site [43].
Acknowledgment.
This work was supported by the projects GA17-12925Sand GA16-03823S of the Czech Science Foundation and by the project LO1506of the Czech Ministry of Education, Youth and Sports.[1] Allaire G (1992) Homogenization and two-scale convergence. SIAM JMath Anal 23:1482–1518, DOI 10.1137/0523084[2] Alnaes MS, Blechta J, Hake J, Johansson A, Kehlet B, Logg A, Richard-son C, Ring J, Rognes ME, Wells GN (2015) The fenics project version1.5. Arch Numer Software 3, DOI 10.11588/ans.2015.100.20553243] Amaya M, Morten JP, Boman L (2016) A low-rank approximation forlarge-scale 3d controlled-source electromagnetic gauss-newton inversion.Geophysics 81(3), DOI 10.1190/geo2015-0079.1[4] Balay S, Abhyankar S, Adams M, Brown J, Brune P, Buschelman K,Dalcin L, Dener A, Eijkhout V, Gropp W, Kaushik D, Knepley M, MayD, Curfman McInnes L, Mills R, Munson T, Rupp K, Sanan P, SmithB, Zampini S, Zhang H, Zhang H (2018) PETSc users manual. Tech.Rep. ANL-95/11 - Revision 3.10, Argonne National Laboratory, URL , accessed 25 September 2018[5] Borden MJ, Scott MA, Evans JA, Hughes TJR (2011) Isogeometricfinite element data structures based on Bezier extraction of NURBS.Int J Numer Meth Engng 87:15–47, DOI 10.1002/nme.2968[6] Bradshaw R, Behnel S, Seljebotn DS, Ewing G, et al. (2018) The Cythoncompiler. Http://cython.org, Accessed 25 September 2018[7] Cimrman R (2014) Enhancing sfepy with isogeometric analysis.arXiv:14126407 [cs, math] URL http://arxiv.org/abs/1412.6407 \\