Initialization

Let us begin by loading Getfem and pyvista.

import getfem as gf
import pyvista as pv
import numpy as np

from pyvirtualdisplay import Display
import os
display = Display(visible=0, size=(1280, 1024))
display.start()

pv.set_plot_theme("document")

Mesh generation

When you want to build mesh from empty situation, you can use “empty” command in Mesh object.

gf.Mesh?

Init signature: gf.Mesh(*args)
Docstring:
GeFEM Mesh object

This object is able to store any element in any dimension even if you mix
elements with different dimensions.
Init docstring:
General constructor for Mesh objects

* ``M = Mesh('empty', int dim)``
  Create a new empty mesh.

* ``M = Mesh('cartesian', vec X[, vec Y[, vec Z,..]])``
  Build quickly a regular mesh of quadrangles, cubes, etc.

* ``M = Mesh('pyramidal', vec X[, vec Y[, vec Z,..]])``
  Build quickly a regular mesh of pyramids, etc.

* ``M = Mesh('cartesian Q1', vec X, vec Y[, vec Z,..])``
  Build quickly a regular mesh of quadrangles, cubes, etc. with
  Q1 elements.

* ``M = Mesh('triangles grid', vec X, vec Y)``
  Build quickly a regular mesh of triangles.

  This is a very limited and somehow deprecated function (See also
  ``Mesh('ptND')``, ``Mesh('regular simplices')`` and
  ``Mesh('cartesian')``).

* ``M = Mesh('regular simplices', vec X[, vec Y[, vec Z,...]]['degree', int k]['noised'])``
  Mesh a n-dimensionnal parallelepipeded with simplices (triangles,
  tetrahedrons etc) .

  The optional degree may be used to build meshes with non linear
  geometric transformations.

* ``M = Mesh('curved', Mesh m, vec F)``
  Build a curved (n+1)-dimensions mesh from a n-dimensions mesh `m`.

  The points of the new mesh have one additional coordinate, given by
  the vector `F`. This can be used to obtain meshes for shells. `m` may
  be a MeshFem object, in that case its linked mesh will be used.

* ``M = Mesh('prismatic', Mesh m, int nl[, int degree])``
  Extrude a prismatic Mesh `M` from a Mesh `m`.

  In the additional dimension there are `nl` layers of elements
  distributed from ``0`` to ``1``.
  If the optional parameter `degree` is provided with a value greater
  than the default value of ``1``, a non-linear transformation of
  corresponding degree is considered in the extrusion direction.

* ``M = Mesh('pt2D', mat P, imat T[, int n])``
  Build a mesh from a 2D triangulation.

  Each column of `P` contains a point coordinate, and each column of `T`
  contains the point indices of a triangle. `n` is optional and is a
  zone number. If `n` is specified then only the zone number `n` is
  converted (in that case, `T` is expected to have 4 rows, the fourth
  containing these zone numbers).



* ``M = Mesh('ptND', mat P, imat T)``
  Build a mesh from a n-dimensional "triangulation".

  Similar function to 'pt2D', for building simplexes meshes from a
  triangulation given in `T`, and a list of points given in `P`. The
  dimension of the mesh will be the number of rows of `P`, and the
  dimension of the simplexes will be the number of rows of `T`.

* ``M = Mesh('load', string filename)``
  Load a mesh from a getfem++ ascii mesh file.

  See also ``Mesh.save(string filename)``.

* ``M = Mesh('from string', string s)``
  Load a mesh from a string description.

  For example, a string returned by ``Mesh.char()``.

* ``M = Mesh('import', string format, string filename)``
  Import a mesh.

  `format` may be:

  - 'gmsh' for a mesh created with `Gmsh`
  - 'gid' for a mesh created with `GiD`
  - 'cdb' for a mesh created with `ANSYS`
  - 'am_fmt' for a mesh created with `EMC2`

* ``M = Mesh('clone', Mesh m2)``
  Create a copy of a mesh.

* ``M = Mesh('generate', MesherObject mo, scalar h[, int K = 1[, mat vertices]])``
  Call the experimental mesher of Getfem on the geometry
  represented by `mo`. please control the conformity of the produced mesh.
  You can help the mesher by adding a priori vertices in the array
  `vertices` which should be of size ``n x m`` where ``n`` n is the
  dimension of the mesh and ``m`` the number of points. `h` is
  approximate diameter of the elements. `K` is the degree of the
  mesh ( > 1 for curved boundaries).  The mesher try to optimize the
  quality of the elements. This operation may be time consuming.
  Note that if the mesh generation fails, because of some random
  procedure used, it can be run again since it will not give necessarily
  the same result due to random procedures used.
  The messages send to the console by the mesh generation can be
  desactivated using `gf_util('trace level', 2)`. More information
  can be obtained by `gf_util('trace level', 4)`. See ``MesherObject``
  to manipulate geometric primitives in order to desribe the geometry.



File:           /usr/lib/python3/dist-packages/getfem/getfem.py
Type:           type
Subclasses:

mesh = gf.Mesh("empty", 3)

for x in range(10):
    for y in range(10):
        for z in [0, 10]:
            mesh.add_convex(
                gf.GeoTrans("GT_PRODUCT(GT_QK(1, 1), GT_QK(1, 1))"),
                [[x, x+1, x, x+1], [y, y, y+1, y+1], [z, z, z, z]],
            )

for x in range(10):
    for z in range(10):
        for y in [0, 10]:
            mesh.add_convex(
                gf.GeoTrans("GT_PRODUCT(GT_QK(1, 1), GT_QK(1, 1))"),
                [[x, x+1, x, x+1], [y, y, y, y], [z, z, z+1, z+1]],
            )

for y in range(10):
    for z in range(10):
        for x in [0, 10]:
            mesh.add_convex(
                gf.GeoTrans("GT_PRODUCT(GT_QK(1, 1), GT_QK(1, 1))"),
                [[x, x, x, x], [y, y+1, y, y+1], [z, z, z+1, z+1]],
           )

mesh.translate([-5.0, -5.0, -5.0])
mesh.export_to_vtk("m1.vtk")

pts = mesh.pts()
radius = np.sqrt(pts[0]**2 + pts[1]**2 + pts[2]**2)
mesh.set_pts(pts*(100.0/radius))
mesh.export_to_vtk("m2.vtk")

p = pv.Plotter(shape=(1, 2))
p.subplot(0, 0)
mesh1 = pv.read("m1.vtk")
p.add_text("Structured mesh", font_size=24)
p.add_mesh(mesh1, show_edges = True)
p.subplot(0, 1)
mesh2 = pv.read("m2.vtk")
p.add_text("Unstructured mesh", font_size=24)
p.add_mesh(mesh2, show_edges = True)
p.show()

[(568.669124079371, 568.669124079371, 568.669124079371),
 (0.0, 0.0, 0.0),
 (0.0, 0.0, 1.0)]

gf.MeshFem?

Init signature: gf.MeshFem(*args)
Docstring:
GeFEM MeshFem object

This object represents a finite element method defined on a whole mesh.
Init docstring:
General constructor for MeshFem objects

* ``MF = MeshFem(Mesh m[, int Qdim1=1[, int Qdim2=1, ...]])``
  Build a new MeshFem object.

  The `Qdim` parameters specifies the dimension of the field represented
  by the finite element method. Qdim1 = 1 for a scalar field,
  Qdim1 = n for a vector field off size n, Qdim1=m, Qdim2=n for
  a matrix field of size mxn ...
  Returns the handle of the created object.

* ``MF = MeshFem('load', string fname[, Mesh m])``
  Load a MeshFem from a file.

  If the mesh `m` is not supplied (this kind of file does not store the
  mesh), then it is read from the file `fname` and its descriptor is
  returned as the second output argument.

* ``MF = MeshFem('from string', string s[, Mesh m])``
  Create a MeshFem object from its string description.

  See also ``MeshFem.char()``

* ``MF = MeshFem('clone', MeshFem mf)``
  Create a copy of a MeshFem.

* ``MF = MeshFem('sum', MeshFem mf1, MeshFem mf2[, MeshFem mf3[, ...]])``
  Create a MeshFem that spans two (or more) MeshFem's.

  All MeshFem must share the same mesh.

  After that, you should not modify the FEM of `mf1`, `mf2` etc.

* ``MF = MeshFem('product', MeshFem mf1, MeshFem mf2)``
  Create a MeshFem that spans all the product of a selection of shape
  functions of `mf1` by all shape functions of `mf2`.
  Designed for Xfem enrichment.

  `mf1` and `mf2` must share the same mesh.

  After that, you should not modify the FEM of `mf1`, `mf2`.

* ``MF = MeshFem('levelset', MeshLevelSet mls, MeshFem mf)``
  Create a MeshFem that is conformal to implicit surfaces defined in
  MeshLevelSet.

* ``MF = MeshFem('global function', Mesh m, LevelSet ls, (GlobalFunction GF1,...)[, int Qdim_m])``
  Create a MeshFem whose base functions are global function given by the
  user in the system of coordinate defined by the iso-values of the two
  level-set function of `ls`.

* ``MF = MeshFem('partial', MeshFem mf, ivec DOFs[, ivec RCVs])``
  Build a restricted MeshFem by keeping only a subset of the degrees of
  freedom of `mf`.

  If `RCVs` is given, no FEM will be put on the convexes listed in
  `RCVs`.


File:           /usr/lib/python3/dist-packages/getfem/getfem.py
Type:           type
Subclasses:

mfu = gf.MeshFem(mesh, 3)
mfd = gf.MeshFem(mesh, 1)
mfu.set_classical_fem(1)
mfd.set_classical_fem(1)
mim = gf.MeshIm(mesh, 2)

model = gf.Model("real")
model.add_fem_variable("u", mfd)
model.add_initialized_data("E", 1.0)
model.add_initialized_data("Nu", 0.3)
model.add_Kirchhoff_Love_plate_brick(mim, "u", "E", "Nu")
model.assembly()
K = model.tangent_matrix()
del model

model = gf.Model("real")
model.add_fem_variable("u", mfd)
model.add_initialized_data("rho", 1.0/(10^6))
model.add_mass_brick(mim, "u", "rho")
model.assembly()
M = model.tangent_matrix()
del model

K.save?

Signature: K.save(format, filename)
Docstring:
Export the sparse matrix.

the format of the file may be 'hb' for Harwell-Boeing, or 'mm'
for Matrix-Market.
File:      /usr/lib/python3/dist-packages/getfem/getfem.py
Type:      method

K.save("mm", "K.mtx")
M.save("mm", "M.mtx")

linalg.eig?

Object `linalg.eig` not found.

import scipy
from scipy import io
from scipy import linalg

SM = io.mmread("K.mtx")
MM = io.mmread("M.mtx")
w, vr = linalg.eig(SM.todense(), MM.todense())
omega = np.sqrt(w.real)

<ipython-input-13-18e62582bc6f>:8: RuntimeWarning: invalid value encountered in sqrt
  omega = np.sqrt(w.real)

p = pv.Plotter(shape=(1, 1))
p.subplot(0, 0)
mesh2 = pv.read("m2.vtk")
p.add_text("Unstructured mesh", font_size=24)
p.add_mesh(mesh2, show_edges = True)
p.show()

[(386.37033051562736, 386.37033051562736, 386.37033051562736),
 (0.0, 0.0, 0.0),
 (0.0, 0.0, 1.0)]

mesh2 = mesh
pts = mesh2.pts()
radius = np.sqrt(pts[0]**2 + pts[1]**2 + pts[2]**2)
mesh2.set_pts(pts*(np.abs(vr[:, 10])/np.max(np.abs(vr[:, 0]))/radius))
mesh2.export_to_vtk("m3.vtk")

p = pv.Plotter(shape=(1, 2))
p.subplot(0, 0)
mesh2 = pv.read("m2.vtk")
p.add_text("Unstructured mesh", font_size=24)
p.add_mesh(mesh2, show_edges = True)
p.subplot(0, 1)
mesh2 = pv.read("m3.vtk")
p.add_text("Eigen mode", font_size=24)
p.add_mesh(mesh2, show_edges = True)
p.show()

[(2.32100872257932, 2.32100872257932, 2.32100872257932),
 (0.0, 0.0, 0.0),
 (0.0, 0.0, 1.0)]