diff --git a/docs/examples/CMakeLists.txt b/docs/examples/CMakeLists.txt
index cb06ae1..9525068 100644
--- a/docs/examples/CMakeLists.txt
+++ b/docs/examples/CMakeLists.txt
@@ -1,60 +1,77 @@
cmake_minimum_required(VERSION 3.0)
if(CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
project(GooseFEM-examples)
find_package(GooseFEM REQUIRED CONFIG)
endif()
option(WARNINGS "Show build warnings" ON)
option(SIMD "Enable xsimd" ON)
option(DEBUG "Enable all assertions" ON)
find_package(HDF5 REQUIRED)
add_library(libraries INTERFACE IMPORTED)
target_link_libraries(libraries INTERFACE GooseFEM ${HDF5_LIBRARIES})
if (SIMD)
target_link_libraries(libraries INTERFACE xtensor::optimize xtensor::use_xsimd)
endif()
if (WARNINGS)
target_link_libraries(libraries INTERFACE GooseFEM::compiler_warnings)
endif()
if (DEBUG)
target_link_libraries(libraries INTERFACE GooseFEM::debug)
endif()
# create executable
set(exec "statics_FixedDisplacements_LinearElastic_example")
set(source "statics/FixedDisplacements_LinearElastic/example.cpp")
add_executable(${exec} ${source})
target_link_libraries(${exec} PRIVATE libraries)
configure_file(
"statics/FixedDisplacements_LinearElastic/example.py"
"statics_FixedDisplacements_LinearElastic_example.py" COPYONLY)
configure_file(
"statics/FixedDisplacements_LinearElastic/plot.py"
"statics_FixedDisplacements_LinearElastic_plot.py" COPYONLY)
# create executable
set(exec "statics_FixedDisplacements_LinearElastic_manual_partition")
set(source "statics/FixedDisplacements_LinearElastic/manual_partition.cpp")
add_executable(${exec} ${source})
target_link_libraries(${exec} PRIVATE libraries)
configure_file(
"statics/FixedDisplacements_LinearElastic/manual_partition.py"
"statics_FixedDisplacements_LinearElastic_manual_partition.py" COPYONLY)
+
+# create executable
+
+set(exec "statics_MixedPeriodic_LinearElastic_example")
+set(source "statics/MixedPeriodic_LinearElastic/example.cpp")
+
+add_executable(${exec} ${source})
+
+target_link_libraries(${exec} PRIVATE libraries)
+
+configure_file(
+ "statics/MixedPeriodic_LinearElastic/example.py"
+ "statics_MixedPeriodic_LinearElastic_example.py" COPYONLY)
+
+configure_file(
+ "statics/MixedPeriodic_LinearElastic/plot.py"
+ "statics_MixedPeriodic_LinearElastic_plot.py" COPYONLY)
diff --git a/docs/examples/statics/MixedPeriodic_LinearElastic/CMakeLists.txt b/docs/examples/statics/MixedPeriodic_LinearElastic/CMakeLists.txt
index 24c53c3..6019030 100644
--- a/docs/examples/statics/MixedPeriodic_LinearElastic/CMakeLists.txt
+++ b/docs/examples/statics/MixedPeriodic_LinearElastic/CMakeLists.txt
@@ -1,64 +1,14 @@
cmake_minimum_required(VERSION 3.1)
-# basic info
-# ----------
-
-# define a project name
-project(main)
-
-# define empty list of libraries to link
-set(PROJECT_LIBS "")
-
-# copy files to build folder
-configure_file("plot.py" "plot.py" COPYONLY)
-
-# basic compiler options
-# ----------------------
-
-# set optimization level
-set(CMAKE_BUILD_TYPE Release)
-
-# set C++ standard
-set(CMAKE_CXX_STANDARD 14)
-set(CMAKE_CXX_STANDARD_REQUIRED ON)
-
-# set warnings on
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -pedantic -g")
-
-# enable optimisation
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native")
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DXTENSOR_USE_XSIMD=ON")
-
-# assertions
-if(NOT ASSERT)
- add_definitions(-DNDEBUG)
-else()
- add_definitions(-DXTENSOR_ENABLE_ASSERT=ON)
-endif()
-
-# load packages
-# -------------
-
-# add custom include paths
-if(NOT "$ENV{INCLUDE_PATH}" STREQUAL "")
- string(REPLACE ":" ";" INCLUDE_LIST "$ENV{INCLUDE_PATH}")
- include_directories(${INCLUDE_LIST})
-endif()
-
-# load HDF5
-find_package(HDF5 COMPONENTS CXX REQUIRED)
-include_directories(${HDF5_INCLUDE_DIRS})
-set(PROJECT_LIBS ${PROJECT_LIBS} ${HDF5_LIBS} ${HDF5_LIBRARIES})
-
-# load eigen
-find_package(PkgConfig)
-pkg_check_modules(EIGEN3 REQUIRED eigen3)
-include_directories(SYSTEM ${EIGEN3_INCLUDE_DIRS})
-
-# create executable
-# -----------------
-
-add_executable(${PROJECT_NAME} main.cpp)
-
-target_link_libraries(${PROJECT_NAME} ${PROJECT_LIBS})
-
+project(example)
+
+set(HIGHFIVE_USE_XTENSOR 1)
+set(HIGHFIVE_USE_BOOST 0)
+find_package(HighFive REQUIRED)
+find_package(xtensor REQUIRED)
+find_package(GooseFEM REQUIRED)
+find_package(GMatElastic REQUIRED)
+
+add_executable(example example.cpp)
+target_link_libraries(example PRIVATE GMatElastic GooseFEM HighFive)
+target_link_libraries(example INTERFACE xtensor::optimize xtensor::use_xsimd)
diff --git a/docs/examples/statics/MixedPeriodic_LinearElastic/example.cpp b/docs/examples/statics/MixedPeriodic_LinearElastic/example.cpp
new file mode 100644
index 0000000..1ac17c7
--- /dev/null
+++ b/docs/examples/statics/MixedPeriodic_LinearElastic/example.cpp
@@ -0,0 +1,145 @@
+#include <GooseFEM/GooseFEM.h>
+#include <GooseFEM/MatrixPartitioned.h>
+#include <GMatElastic/Cartesian3d.h>
+#include <highfive/H5Easy.hpp>
+
+int main()
+{
+ // mesh
+ // ----
+
+ // define mesh
+ GooseFEM::Mesh::Quad4::Regular mesh(5, 5);
+
+ // mesh dimensions
+ size_t nelem = mesh.nelem();
+ size_t nne = mesh.nne();
+ size_t ndim = mesh.ndim();
+
+ // mesh definitions
+ xt::xtensor<double,2> coor = mesh.coor();
+ xt::xtensor<size_t,2> conn = mesh.conn();
+ xt::xtensor<size_t,2> dofs = mesh.dofs();
+
+ // node sets
+ xt::xtensor<size_t,1> nodesLft = mesh.nodesLeftOpenEdge();
+ xt::xtensor<size_t,1> nodesRgt = mesh.nodesRightOpenEdge();
+ xt::xtensor<size_t,1> nodesTop = mesh.nodesTopEdge();
+ xt::xtensor<size_t,1> nodesBot = mesh.nodesBottomEdge();
+
+ // periodicity and fixed displacements DOFs
+ // ----------------------------------------
+
+ for (size_t j = 0; j < coor.shape(1); ++j) {
+ xt::view(dofs, xt::keep(nodesRgt), j) = xt::view(dofs, xt::keep(nodesLft), j);
+ }
+
+ dofs = GooseFEM::Mesh::renumber(dofs);
+
+ xt::xtensor<size_t,1> iip = xt::concatenate(xt::xtuple(
+ xt::view(dofs, xt::keep(nodesBot), 0),
+ xt::view(dofs, xt::keep(nodesBot), 1),
+ xt::view(dofs, xt::keep(nodesTop), 0),
+ xt::view(dofs, xt::keep(nodesTop), 1)));
+
+ // simulation variables
+ // --------------------
+
+ // vector definition
+ GooseFEM::VectorPartitioned vector(conn, dofs, iip);
+
+ // allocate system matrix
+ GooseFEM::MatrixPartitioned<> K(conn, dofs, iip);
+
+ // nodal quantities
+ xt::xtensor<double,2> disp = xt::zeros<double>(coor.shape());
+ xt::xtensor<double,2> fint = xt::zeros<double>(coor.shape());
+ xt::xtensor<double,2> fext = xt::zeros<double>(coor.shape());
+ xt::xtensor<double,2> fres = xt::zeros<double>(coor.shape());
+
+ // element vectors
+ xt::xtensor<double,3> ue = xt::empty<double>({nelem, nne, ndim});
+ xt::xtensor<double,3> fe = xt::empty<double>({nelem, nne, ndim});
+ xt::xtensor<double,3> Ke = xt::empty<double>({nelem, nne * ndim, nne * ndim});
+
+ // element/material definition
+ // ---------------------------
+
+ // element definition
+ GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor));
+ size_t nip = elem.nip();
+
+ // material definition
+ GMatElastic::Cartesian3d::Matrix mat(nelem, nip);
+ xt::xtensor<size_t,2> Ihard = xt::zeros<size_t>({nelem, nip});
+ xt::view(Ihard, xt::keep(0, 1, 5, 6), xt::all()) = 1;
+ xt::xtensor<size_t,2> Isoft = xt::ones<size_t>({nelem, nip}) - Ihard;
+ mat.setElastic(Isoft, 10.0, 1.0);
+ mat.setElastic(Ihard, 10.0, 10.0);
+
+ // integration point tensors
+ xt::xtensor<double,4> Eps = xt::empty<double>({nelem, nip, 3ul, 3ul});
+ xt::xtensor<double,4> Sig = xt::empty<double>({nelem, nip, 3ul, 3ul});
+ xt::xtensor<double,6> C = xt::empty<double>({nelem, nip, 3ul, 3ul, 3ul, 3ul});
+
+ // solve
+ // -----
+
+ // strain
+ vector.asElement(disp, ue);
+ elem.symGradN_vector(ue, Eps);
+
+ // stress & tangent
+ mat.tangent(Eps, Sig, C);
+
+ // internal force
+ elem.int_gradN_dot_tensor2_dV(Sig, fe);
+ vector.assembleNode(fe, fint);
+
+ // stiffness matrix
+ elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke);
+ K.assemble(Ke);
+
+ // set fixed displacements
+ xt::view(disp, xt::keep(nodesTop), 0) = +0.1;
+
+ // residual
+ xt::noalias(fres) = fext - fint;
+
+ // solve
+ K.solve(fres, disp);
+
+ // post-process
+ // ------------
+
+ // compute strain and stress
+ vector.asElement(disp, ue);
+ elem.symGradN_vector(ue, Eps);
+ mat.stress(Eps, Sig);
+
+ // internal force
+ elem.int_gradN_dot_tensor2_dV(Sig, fe);
+ vector.assembleNode(fe, fint);
+
+ // apply reaction force
+ vector.copy_p(fint, fext);
+
+ // residual
+ xt::noalias(fres) = fext - fint;
+
+ // print residual
+ std::cout << xt::sum(xt::abs(fres))[0] / xt::sum(xt::abs(fext))[0] << std::endl;
+
+ // average stress per node
+ xt::xtensor<double,4> dV = elem.DV(2);
+ xt::xtensor<double,3> SigAv = xt::average(Sig, dV, {1});
+
+ // write output
+ H5Easy::File file("output.h5", H5Easy::File::Overwrite);
+ H5Easy::dump(file, "/coor", coor);
+ H5Easy::dump(file, "/conn", conn);
+ H5Easy::dump(file, "/disp", disp);
+ H5Easy::dump(file, "/Sig", SigAv);
+
+ return 0;
+}
diff --git a/docs/examples/statics/MixedPeriodic_LinearElastic/example.py b/docs/examples/statics/MixedPeriodic_LinearElastic/example.py
new file mode 100644
index 0000000..d56432a
--- /dev/null
+++ b/docs/examples/statics/MixedPeriodic_LinearElastic/example.py
@@ -0,0 +1,188 @@
+import sys
+import GooseFEM
+import GMatElastic
+import numpy as np
+
+# mesh
+# ----
+
+# define mesh
+mesh = GooseFEM.Mesh.Quad4.Regular(5, 5)
+
+# mesh dimensions
+nelem = mesh.nelem()
+nne = mesh.nne()
+ndim = mesh.ndim()
+
+# mesh definitions
+coor = mesh.coor()
+conn = mesh.conn()
+dofs = mesh.dofs()
+
+# node sets
+nodesLft = mesh.nodesLeftOpenEdge()
+nodesRgt = mesh.nodesRightOpenEdge()
+nodesTop = mesh.nodesTopEdge()
+nodesBot = mesh.nodesBottomEdge()
+
+# periodicity and fixed displacements DOFs
+# ----------------------------------------
+
+for j in range(coor.shape[1]):
+ dofs[nodesRgt, j] = dofs[nodesLft, j]
+
+dofs = GooseFEM.Mesh.renumber(dofs)
+
+iip = np.concatenate((
+ dofs[nodesBot, 0],
+ dofs[nodesBot, 1],
+ dofs[nodesTop, 0],
+ dofs[nodesTop, 1]
+))
+
+# simulation variables
+# --------------------
+
+# vector definition
+vector = GooseFEM.VectorPartitioned(conn, dofs, iip)
+
+# allocate system matrix
+K = GooseFEM.MatrixPartitioned(conn, dofs, iip)
+
+# nodal quantities
+disp = np.zeros(coor.shape)
+fint = np.zeros(coor.shape)
+fext = np.zeros(coor.shape)
+fres = np.zeros(coor.shape)
+
+# element vectors
+ue = np.empty((nelem, nne, ndim))
+fe = np.empty((nelem, nne, ndim))
+Ke = np.empty((nelem, nne * ndim, nne * ndim))
+
+# element/material definition
+# ---------------------------
+
+# element definition
+elem = GooseFEM.Element.Quad4.QuadraturePlanar(vector.AsElement(coor))
+nip = elem.nip()
+
+# material definition
+mat = GMatElastic.Cartesian3d.Matrix(nelem, nip)
+Ihard = np.zeros((nelem, nip), dtype='int')
+Ihard[[0, 1, 5, 6], :] = 1
+Isoft = np.ones((nelem, nip), dtype='int') - Ihard
+mat.setElastic(Isoft, 10.0, 1.0)
+mat.setElastic(Ihard, 10.0, 10.0)
+
+# integration point tensors
+Eps = np.empty((nelem, nip, 3, 3))
+Sig = np.empty((nelem, nip, 3, 3))
+C = np.empty((nelem, nip, 3, 3, 3, 3))
+
+# solve
+# -----
+
+# strain
+ue = vector.AsElement(disp)
+Eps = elem.SymGradN_vector(ue)
+
+# stress & tangent
+(Sig, C) = mat.Tangent(Eps)
+
+# internal force
+fe = elem.Int_gradN_dot_tensor2_dV(Sig)
+fint = vector.AssembleNode(fe)
+
+# stiffness matrix
+Ke = elem.Int_gradN_dot_tensor4_dot_gradNT_dV(C)
+K.assemble(Ke)
+
+# set fixed displacements
+disp[nodesTop, 0] = +0.1
+
+# residual
+fres = fext - fint
+
+# solve
+disp = K.Solve(fres, disp)
+
+# post-process
+# ------------
+
+# compute strain and stress
+ue = vector.AsElement(disp)
+Eps = elem.SymGradN_vector(ue)
+Sig = mat.Stress(Eps)
+
+# internal force
+fe = elem.Int_gradN_dot_tensor2_dV(Sig)
+fint = vector.AssembleNode(fe)
+
+# apply reaction force
+fext = vector.Copy_p(fint, fext)
+
+# residual
+fres = fext - fint
+
+# print residual
+print(np.sum(np.abs(fres)) / np.sum(np.abs(fext)))
+
+# average stress per node
+dV = elem.DV(2)
+Sig = np.average(Sig, weights=dV, axis=1)
+
+# skip plot with "--no-plot" command line argument
+# ------------------------------------------------
+
+if len(sys.argv) == 2:
+ if sys.argv[1] == "--no-plot":
+ sys.exit(0)
+
+# plot
+# ----
+
+import matplotlib.pyplot as plt
+import GooseMPL as gplt
+
+plt.style.use(['goose', 'goose-latex'])
+
+# extract dimension
+
+nelem = conn.shape[0]
+
+# tensor products
+
+def ddot22(A2, B2):
+ return np.einsum('eij, eji -> e', A2, B2)
+
+
+def ddot42(A4, B2):
+ return np.einsum('eijkl, elk -> eij', A4, B2)
+
+
+def dyad22(A2, B2):
+ return np.einsum('eij, ekl -> eijkl', A2, B2)
+
+# identity tensors
+
+i = np.eye(3)
+I = np.einsum('ij, e', i, np.ones([nelem]))
+I4 = np.einsum('ijkl, e -> eijkl', np.einsum('il, jk', i, i), np.ones([nelem]))
+I4rt = np.einsum('ijkl, e -> eijkl', np.einsum('ik,jl', i, i), np.ones([nelem]))
+I4s = 0.5 * (I4 + I4rt)
+II = dyad22(I, I)
+I4d = I4s - II / 3.0
+
+# compute equivalent stress
+
+Sigd = ddot42(I4d, Sig)
+sigeq = np.sqrt(3.0 / 2.0 * ddot22(Sigd, Sigd))
+
+# plot
+
+fig, ax = plt.subplots()
+gplt.patch(coor=coor + disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0, 0.1))
+gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax)
+plt.savefig('plot.pdf')
+plt.close()
diff --git a/docs/examples/statics/MixedPeriodic_LinearElastic/main.cpp b/docs/examples/statics/MixedPeriodic_LinearElastic/main.cpp
deleted file mode 100644
index 8f53f81..0000000
--- a/docs/examples/statics/MixedPeriodic_LinearElastic/main.cpp
+++ /dev/null
@@ -1,159 +0,0 @@
-#include <Eigen/Eigen>
-#include <GooseFEM/GooseFEM.h>
-#include <GMatElastic/Cartesian3d.h>
-#include <xtensor-io/xhighfive.hpp>
-
-int main()
-{
- // mesh
- // ----
-
- // define mesh
-
- GooseFEM::Mesh::Quad4::Regular mesh(5,5);
-
- // mesh dimensions
-
- size_t nelem = mesh.nelem();
- size_t nne = mesh.nne();
- size_t ndim = mesh.ndim();
-
- // mesh definitions
-
- xt::xtensor<double,2> coor = mesh.coor();
- xt::xtensor<size_t,2> conn = mesh.conn();
- xt::xtensor<size_t,2> dofs = mesh.dofs();
-
- // node sets
-
- xt::xtensor<size_t,1> nodesLeft = mesh.nodesLeftOpenEdge();
- xt::xtensor<size_t,1> nodesRight = mesh.nodesRightOpenEdge();
- xt::xtensor<size_t,1> nodesTop = mesh.nodesTopEdge();
- xt::xtensor<size_t,1> nodesBottom = mesh.nodesBottomEdge();
-
- // periodicity and fixed displacements DOFs
- // ----------------------------------------
-
- for ( size_t i = 0 ; i < nodesLeft.size() ; ++i )
- for ( size_t j = 0 ; j < coor.shape()[1] ; ++j )
- dofs(nodesRight(i),j) = dofs(nodesLeft(i),j);
-
- dofs = GooseFEM::Mesh::renumber(dofs);
-
- xt::xtensor<size_t,1> iip = xt::concatenate(xt::xtuple(
- xt::view(dofs, xt::keep(nodesBottom), 0),
- xt::view(dofs, xt::keep(nodesBottom), 1),
- xt::view(dofs, xt::keep(nodesTop ), 0),
- xt::view(dofs, xt::keep(nodesTop ), 1)
- ));
-
- // simulation variables
- // --------------------
-
- // vector definition
-
- GooseFEM::VectorPartitioned vector(conn, dofs, iip);
-
- // nodal quantities
-
- xt::xtensor<double,2> disp = xt::zeros<double>(coor.shape());
- xt::xtensor<double,2> fint = xt::zeros<double>(coor.shape());
- xt::xtensor<double,2> fext = xt::zeros<double>(coor.shape());
- xt::xtensor<double,2> fres = xt::zeros<double>(coor.shape());
-
- // element vectors
-
- xt::xtensor<double,3> ue = xt::empty<double>({nelem, nne, ndim});
- xt::xtensor<double,3> fe = xt::empty<double>({nelem, nne, ndim});
- xt::xtensor<double,3> Ke = xt::empty<double>({nelem, nne*ndim, nne*ndim});
-
- // element/material definition
- // ---------------------------
-
- GooseFEM::Element::Quad4::QuadraturePlanar elem(vector.AsElement(coor));
-
- size_t nip = elem.nip();
-
- GMatElastic::Cartesian3d::Matrix mat(nelem, nip);
-
- xt::xtensor<size_t,2> Ihard = xt::zeros<size_t>({nelem, nip});
-
- xt::view(Ihard, xt::keep(0,1,5,6), xt::all()) = 1;
-
- xt::xtensor<size_t,2> Isoft = xt::ones<size_t>({nelem, nip}) - Ihard;
-
- mat.setElastic(Isoft,10.,1. );
- mat.setElastic(Ihard,10.,10.);
-
- // solve
- // -----
-
- // allocate tensors
- size_t d = 3;
- xt::xtensor<double,4> Eps = xt::empty<double>({nelem, nip, d, d });
- xt::xtensor<double,4> Sig = xt::empty<double>({nelem, nip, d, d });
- xt::xtensor<double,6> C = xt::empty<double>({nelem, nip, d, d, d, d});
-
- // allocate system matrix
- GooseFEM::MatrixPartitioned<> K(conn, dofs, iip);
-
- // strain
- vector.asElement(disp, ue);
- elem.symGradN_vector(ue, Eps);
-
- // stress & tangent
- mat.tangent(Eps, Sig, C);
-
- // internal force
- elem.int_gradN_dot_tensor2_dV(Sig, fe);
- vector.assembleNode(fe, fint);
-
- // stiffness matrix
- elem.int_gradN_dot_tensor4_dot_gradNT_dV(C, Ke);
- K.assemble(Ke);
-
- // set fixed displacements
- xt::view(disp, xt::keep(nodesTop), xt::keep(0)) = +0.1;
-
- // residual
- xt::noalias(fres) = fext - fint;
-
- // solve
- K.solve(fres, disp);
-
- // post-process
- // ------------
-
- // compute strain and stress
- vector.asElement(disp, ue);
- elem.symGradN_vector(ue, Eps);
- mat.stress(Eps, Sig);
-
- // internal force
- elem.int_gradN_dot_tensor2_dV(Sig, fe);
- vector.assembleNode(fe, fint);
-
- // apply reaction force
- vector.copy_p(fint, fext);
-
- // residual
- xt::noalias(fres) = fext - fint;
-
- // print residual
- std::cout << xt::sum(xt::abs(fres))[0] / xt::sum(xt::abs(fext))[0] << std::endl;
-
- // average stress per node
- xt::xtensor<double,4> dV = elem.DV(2);
- xt::xtensor<double,3> SigAv = xt::average(Sig, dV, {1});
-
- // write output
-
- HighFive::File file("main.h5", HighFive::File::Overwrite);
-
- xt::dump(file, "/coor", coor);
- xt::dump(file, "/conn", conn);
- xt::dump(file, "/disp", disp);
- xt::dump(file, "/Sig" , SigAv);
-
- return 0;
-}
diff --git a/docs/examples/statics/MixedPeriodic_LinearElastic/plot.py b/docs/examples/statics/MixedPeriodic_LinearElastic/plot.py
index 0bf1209..8447c77 100644
--- a/docs/examples/statics/MixedPeriodic_LinearElastic/plot.py
+++ b/docs/examples/statics/MixedPeriodic_LinearElastic/plot.py
@@ -1,51 +1,51 @@
-
import h5py
-
import matplotlib.pyplot as plt
-import GooseMPL as gplt
-import numpy as np
+import GooseMPL as gplt
+import numpy as np
plt.style.use(['goose', 'goose-latex'])
-# open file
-file = h5py.File('main.h5', 'r')
+# load data
-# read fields
-coor = file['/coor'][...]
-conn = file['/conn'][...]
-disp = file['/disp'][...]
-Sig = file['/Sig' ][...]
-
-# extract dimension
-nelem = conn.shape[0]
+with h5py.File('output.h5', 'r') as data:
+ coor = data['/coor'][...]
+ conn = data['/conn'][...]
+ disp = data['/disp'][...]
+ Sig = data['/Sig'][...]
+ nelem = conn.shape[0]
# tensor products
-ddot22 = lambda A2,B2: np.einsum('eij ,eji->e ',A2,B2)
-ddot42 = lambda A4,B2: np.einsum('eijkl,elk->eij ',A4,B2)
-dyad22 = lambda A2,B2: np.einsum('eij ,ekl->eijkl',A2,B2)
-# identity tensor (single tensor)
-i = np.eye(3)
+def ddot22(A2, B2):
+ return np.einsum('eij, eji -> e', A2, B2)
-# identity tensors (grid)
-I = np.einsum('ij ,e' , i ,np.ones([nelem]))
-I4 = np.einsum('ijkl,e->eijkl',np.einsum('il,jk',i,i),np.ones([nelem]))
-I4rt = np.einsum('ijkl,e->eijkl',np.einsum('ik,jl',i,i),np.ones([nelem]))
-I4s = (I4+I4rt)/2.
-II = dyad22(I,I)
-I4d = I4s-II/3.
-# compute equivalent stress
-Sigd = ddot42(I4d, Sig)
-sigeq = np.sqrt(3./2.*ddot22(Sigd,Sigd))
+def ddot42(A4, B2):
+ return np.einsum('eijkl, elk -> eij', A4, B2)
-# plot
-fig, ax = plt.subplots()
+def dyad22(A2, B2):
+ return np.einsum('eij, ekl -> eijkl', A2, B2)
-gplt.patch(coor=coor+disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0,0.1))
+# identity tensors
-gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax)
+i = np.eye(3)
+I = np.einsum('ij, e', i, np.ones([nelem]))
+I4 = np.einsum('ijkl, e -> eijkl', np.einsum('il, jk', i, i), np.ones([nelem]))
+I4rt = np.einsum('ijkl, e -> eijkl', np.einsum('ik,jl', i, i), np.ones([nelem]))
+I4s = 0.5 * (I4 + I4rt)
+II = dyad22(I, I)
+I4d = I4s - II / 3.0
-plt.show()
+# compute equivalent stress
+
+Sigd = ddot42(I4d, Sig)
+sigeq = np.sqrt(3.0 / 2.0 * ddot22(Sigd, Sigd))
+
+# plot
+fig, ax = plt.subplots()
+gplt.patch(coor=coor + disp, conn=conn, cindex=sigeq, cmap='jet', axis=ax, clim=(0, 0.1))
+gplt.patch(coor=coor, conn=conn, linestyle='--', axis=ax)
+plt.savefig('plot.pdf')
+plt.close()