finsberg · March 8, 2024 12:45 · finsberg · Mar 8, 2024
diff --git a/Dockerfile_fp32 b/Dockerfile_fp32
 from ghcr.io/fenics/test-env:current-mpich

 ENV PETSC_ARCH=linux-gnu-real32-32
 ARG UFL_BRANCH=main
 ARG FFCX_BRANCH=main
 ARG BUILD_MODE=Release

 RUN python3 -m pip install pybind11

 RUN git clone --branch=main --single-branch https://github.com/FEniCS/dolfinx.git
 RUN git clone --branch=main --single-branch https://github.com/FEniCS/basix.git
 RUN cmake -G Ninja -B build-basix -DCMAKE_BUILD_TYPE="Release" -S ./basix/cpp/
 RUN cmake --build build-basix --parallel 2
 RUN cmake --install build-basix
 RUN python3 -m pip -v install --check-build-dependencies --config-settings=build-dir="build-basix-python" --config-settings=cmake.build-type=${BUILD_MODE} --config-settings=install.strip=false --no-build-isolation ./basix/python
 RUN python3 -m pip install git+https://github.com/FEniCS/ufl.git@${UFL_BRANCH}
 RUN python3 -m pip install git+https://github.com/FEniCS/ffcx.git@${FFCX_BRANCH}


 RUN git clone --branch=master --single-branch https://github.com/ornladios/ADIOS2.git adios2
 RUN cmake -G Ninja -DADIOS2_USE_HDF5=on -DCMAKE_INSTALL_PYTHONDIR=/usr/local/lib/ -DADIOS2_USE_Fortran=off -DBUILD_TESTING=off -DADIOS2_BUILD_EXAMPLES=off -DADIOS2_USE_ZeroMQ=off -B build-dir-adios2 -S ./adios2
 RUN cmake -G Ninja -B build-dir-adios2 -DCMAKE_BUILD_TYPE="Release" -S ./adios2/
 RUN cmake --build build-dir-adios2 --parallel 2
 RUN cmake --install build-dir-adios2

 ENV PETSC_DIR=/usr/local/petsc
 RUN PETSC_DIR=$PETSC_DIR PETSC_ARCH=${PETSC_ARCH} cmake -G Ninja -DCMAKE_BUILD_TYPE=${BUILD_MODE} -B build-dolfinx -S ./dolfinx/cpp/

 RUN cmake --build build-dolfinx
 RUN cmake --install build-dolfinx

 RUN PETSC_DIR=${PETSC_DIR} PETSC_ARCH=${PETSC_ARCH} python3 -m pip -v install --check-build-dependencies --no-build-isolation  ./dolfinx/python/
diff --git a/hyperelasticity.py b/hyperelasticity.py
 from dolfinx import log, default_scalar_type
 from dolfinx.fem.petsc import NonlinearProblem
 from dolfinx.nls.petsc import NewtonSolver
 from dolfinx import nls
 import numpy as np
 import ufl
 import dolfinx

 from mpi4py import MPI
 from dolfinx import fem, mesh, plot
 L = 20.0
 domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron, dtype=dolfinx.default_real_type)
 V = dolfinx.fem.functionspace(domain, ("Lagrange", 2, (3,)))


 # +
 def left(x):
    return np.isclose(x[0], 0)


 def right(x):
    return np.isclose(x[0], L)


 fdim = domain.topology.dim - 1
 left_facets = mesh.locate_entities_boundary(domain, fdim, left)
 right_facets = mesh.locate_entities_boundary(domain, fdim, right)
 # -

 # Next, we create a  marker based on these two functions

 # Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
 marked_facets = np.hstack([left_facets, right_facets])
 marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
 sorted_facets = np.argsort(marked_facets)
 facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

 # We then create a function for supplying the boundary condition on the left side, which is fixed.

 u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)

 # To apply the boundary condition, we identity the dofs located on the facets marked by the `MeshTag`.

 left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
 bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

 # Next, we define the body force on the reference configuration (`B`), and nominal (first Piola-Kirchhoff) traction (`T`).

 B = fem.Constant(domain, default_scalar_type((0, 0, 0)))
 T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

 # Define the test and solution functions on the space $V$

 v = ufl.TestFunction(V)
 u = fem.Function(V)

 # +
 # Spatial dimension
 d = len(u)

 # Identity tensor
 I = ufl.variable(ufl.Identity(d))

 # Deformation gradient
 F = ufl.variable(I + ufl.grad(u))

 # Right Cauchy-Green tensor
 C = ufl.variable(F.T * F)

 # Invariants of deformation tensors
 Ic = ufl.variable(ufl.tr(C))
 J = ufl.variable(ufl.det(F))
 # -

 # Define the elasticity model via a stored strain energy density function $\psi$, and create the expression for the first Piola-Kirchhoff stress:

 # Elasticity parameters
 E = default_scalar_type(1.0e4)
 nu = default_scalar_type(0.3)
 mu = fem.Constant(domain, default_scalar_type(E / (2 * (1 + nu))))
 lmbda = fem.Constant(domain, default_scalar_type(E * nu / ((1 + nu) * (1 - 2 * nu))))
 # Stored strain energy density (compressible neo-Hookean model)
 psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
 # Stress
 # Hyper-elasticity
 P = ufl.diff(psi, F)

 metadata = {"quadrature_degree": 6}
 ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
 dx = ufl.Measure("dx", domain=domain)#, metadata=metadata)
 # Define form F (we want to find u such that F(u) = 0)
 F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)

 # As the varitional form is non-linear and written on residual form, we use the non-linear problem class from DOLFINx to set up required structures to use a Newton solver.

 problem = NonlinearProblem(F, u, bcs)

 # and then create and customize the Newton solver

 # +
 solver = NewtonSolver(domain.comm, problem)

 # Set Newton solver options
 solver.atol = 1e-5
 solver.rtol = 1e-5
 solver.krylov_solver.setType("preonly")
 solver.krylov_solver.getPC().setType("lu")
 solver.krylov_solver.getPC().setFactorSolverType("mumps")

 solver.max_it = 25
 solver.convergence_criterion = "residual"


 log.set_log_level(log.LogLevel.INFO)
 tval0 = -1.5
 for n in range(1, 2):
    T.value[2] = n * tval0
    num_its, converged = solver.solve(u)
    assert (converged)
    u.x.scatter_forward()
    print(f"Time step {n}, Number of iterations {num_its}, Load {T.value}")
	from ghcr.io/fenics/test-env:current-mpich

	ENV PETSC_ARCH=linux-gnu-real32-32
	ARG UFL_BRANCH=main
	ARG FFCX_BRANCH=main
	ARG BUILD_MODE=Release

	RUN python3 -m pip install pybind11

	RUN git clone --branch=main --single-branch https://github.com/FEniCS/dolfinx.git
	RUN git clone --branch=main --single-branch https://github.com/FEniCS/basix.git
	RUN cmake -G Ninja -B build-basix -DCMAKE_BUILD_TYPE="Release" -S ./basix/cpp/
	RUN cmake --build build-basix --parallel 2
	RUN cmake --install build-basix
	RUN python3 -m pip -v install --check-build-dependencies --config-settings=build-dir="build-basix-python" --config-settings=cmake.build-type=${BUILD_MODE} --config-settings=install.strip=false --no-build-isolation ./basix/python
	RUN python3 -m pip install git+https://github.com/FEniCS/ufl.git@${UFL_BRANCH}
	RUN python3 -m pip install git+https://github.com/FEniCS/ffcx.git@${FFCX_BRANCH}


	RUN git clone --branch=master --single-branch https://github.com/ornladios/ADIOS2.git adios2
	RUN cmake -G Ninja -DADIOS2_USE_HDF5=on -DCMAKE_INSTALL_PYTHONDIR=/usr/local/lib/ -DADIOS2_USE_Fortran=off -DBUILD_TESTING=off -DADIOS2_BUILD_EXAMPLES=off -DADIOS2_USE_ZeroMQ=off -B build-dir-adios2 -S ./adios2
	RUN cmake -G Ninja -B build-dir-adios2 -DCMAKE_BUILD_TYPE="Release" -S ./adios2/
	RUN cmake --build build-dir-adios2 --parallel 2
	RUN cmake --install build-dir-adios2

	ENV PETSC_DIR=/usr/local/petsc
	RUN PETSC_DIR=$PETSC_DIR PETSC_ARCH=${PETSC_ARCH} cmake -G Ninja -DCMAKE_BUILD_TYPE=${BUILD_MODE} -B build-dolfinx -S ./dolfinx/cpp/

	RUN cmake --build build-dolfinx
	RUN cmake --install build-dolfinx

	RUN PETSC_DIR=${PETSC_DIR} PETSC_ARCH=${PETSC_ARCH} python3 -m pip -v install --check-build-dependencies --no-build-isolation ./dolfinx/python/
	from dolfinx import log, default_scalar_type
	from dolfinx.fem.petsc import NonlinearProblem
	from dolfinx.nls.petsc import NewtonSolver
	from dolfinx import nls
	import numpy as np
	import ufl
	import dolfinx

	from mpi4py import MPI
	from dolfinx import fem, mesh, plot
	L = 20.0
	domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron, dtype=dolfinx.default_real_type)
	V = dolfinx.fem.functionspace(domain, ("Lagrange", 2, (3,)))


	# +
	def left(x):
	return np.isclose(x[0], 0)


	def right(x):
	return np.isclose(x[0], L)


	fdim = domain.topology.dim - 1
	left_facets = mesh.locate_entities_boundary(domain, fdim, left)
	right_facets = mesh.locate_entities_boundary(domain, fdim, right)
	# -

	# Next, we create a marker based on these two functions

	# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
	marked_facets = np.hstack([left_facets, right_facets])
	marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
	sorted_facets = np.argsort(marked_facets)
	facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

	# We then create a function for supplying the boundary condition on the left side, which is fixed.

	u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)

	# To apply the boundary condition, we identity the dofs located on the facets marked by the `MeshTag`.

	left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
	bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

	# Next, we define the body force on the reference configuration (`B`), and nominal (first Piola-Kirchhoff) traction (`T`).

	B = fem.Constant(domain, default_scalar_type((0, 0, 0)))
	T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

	# Define the test and solution functions on the space $V$

	v = ufl.TestFunction(V)
	u = fem.Function(V)

	# +
	# Spatial dimension
	d = len(u)

	# Identity tensor
	I = ufl.variable(ufl.Identity(d))

	# Deformation gradient
	F = ufl.variable(I + ufl.grad(u))

	# Right Cauchy-Green tensor
	C = ufl.variable(F.T * F)

	# Invariants of deformation tensors
	Ic = ufl.variable(ufl.tr(C))
	J = ufl.variable(ufl.det(F))
	# -

	# Define the elasticity model via a stored strain energy density function $\psi$, and create the expression for the first Piola-Kirchhoff stress:

	# Elasticity parameters
	E = default_scalar_type(1.0e4)
	nu = default_scalar_type(0.3)
	mu = fem.Constant(domain, default_scalar_type(E / (2 * (1 + nu))))
	lmbda = fem.Constant(domain, default_scalar_type(E * nu / ((1 + nu) * (1 - 2 * nu))))
	# Stored strain energy density (compressible neo-Hookean model)
	psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
	# Stress
	# Hyper-elasticity
	P = ufl.diff(psi, F)

	metadata = {"quadrature_degree": 6}
	ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
	dx = ufl.Measure("dx", domain=domain)#, metadata=metadata)
	# Define form F (we want to find u such that F(u) = 0)
	F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)

	# As the varitional form is non-linear and written on residual form, we use the non-linear problem class from DOLFINx to set up required structures to use a Newton solver.

	problem = NonlinearProblem(F, u, bcs)

	# and then create and customize the Newton solver

	# +
	solver = NewtonSolver(domain.comm, problem)

	# Set Newton solver options
	solver.atol = 1e-5
	solver.rtol = 1e-5
	solver.krylov_solver.setType("preonly")
	solver.krylov_solver.getPC().setType("lu")
	solver.krylov_solver.getPC().setFactorSolverType("mumps")

	solver.max_it = 25
	solver.convergence_criterion = "residual"


	log.set_log_level(log.LogLevel.INFO)
	tval0 = -1.5
	for n in range(1, 2):
	T.value[2] = n * tval0
	num_its, converged = solver.solve(u)
	assert (converged)
	u.x.scatter_forward()
	print(f"Time step {n}, Number of iterations {num_its}, Load {T.value}")