Set initial conditions for tensor-element

dolfinx
1

Hello, I have a system with mixed elements; 2 vector elements and 2 tensor elements (2-dimensional tensors). Since I have a system, where one of these tensors is invertet in the residuum, I need to start with a non-singular tensor, for example the identity matrix.

from dolfinx import *
from mpi4py import MPI
import numpy as np
import matplotlib.pyplot as plt
import ufl
import scipy.io
from petsc4py import PETSc

mesh = mesh.create_unit_cube(MPI.COMM_WORLD, nx=1, ny=1, nz=1)
n = ufl.FacetNormal(mesh)

P = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
R = ufl.TensorElement("DG", mesh.ufl_cell(), 1)
P_wish = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
R_wish = ufl.TensorElement("DG", mesh.ufl_cell(), 1)
mel = ufl.MixedElement([P, R, P_wish, R_wish])
print(mel)
U = fem.FunctionSpace(mesh, mel)

u = fem.Function(U)
(p, F, v, P) = ufl.split(u)
(dv, dP, dp, dF) = ufl.TestFunctions(U)

rho_0=1
mu=1
lamb=1

Psi_fun=lambda F: mu/2*(np.trace(F.T@F)-3)-mu*np.log(np.sqrt(np.linalg.det(F.T@F)))+lamb/2*np.log(np.sqrt(np.linalg.det(F.T@F)))**2
P_fun=lambda F: ufl.dot(F,mu*(ufl.classes.Identity(F.ufl_shape[0])-ufl.inv(ufl.dot(F.T,F)))+lamb*ufl.ln(ufl.sqrt(ufl.classes.Determinant(ufl.dot(F.T,F))))*ufl.inv(ufl.dot(F.T,F)))

u_n = fem.Function(U) # previous time step
(p_n, F_n, v_n, P_n) = ufl.split(u_n)

# Initial conditions:

def initial_cond(x,y,z):
    return np.eye(F.ufl_shape[0])
F.interpolate(np.eye(F.ufl_shape[0]))
F_n.interpolate(np.eye(F.ufl_shape[0]))

It just tells me

AttributeError: 'ListTensor' object has no attribute 'interpolate'

so what I found on in the examples does not work, unfortunately. Can anybody help me how to set a tensor-function-space initially to identity?

Thank you in advance for any help!

2

To set initial conditions, you should use:
u.split() or u.sub(1).interpolate
while `ufl.split(u) should be used when declaring variables in the linear and/or bilinear form

3

So, even if I already called it “F”, I now add the command

u.sub(1).interpolate(initial_cond)

instead of

F.interpolate(np.eye(F.ufl_shape[0]))

or how do you mean it? I don’t quite have the feeling yet where to apply this. Or do I do this when I first split u?

4

Yes, consider the minimal example:

import dolfinx
from mpi4py import MPI
import numpy as np
import ufl

mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, nx=1, ny=1, nz=1)
n = ufl.FacetNormal(mesh)

P = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
R = ufl.TensorElement("DG", mesh.ufl_cell(), 1)
P_wish = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
R_wish = ufl.TensorElement("DG", mesh.ufl_cell(), 1)
mel = ufl.MixedElement([P, R, P_wish, R_wish])
U = dolfinx.fem.FunctionSpace(mesh, mel)

u = dolfinx.fem.Function(U)
(p, F, v, P) = ufl.split(u)


def T(x):
    values = np.zeros((mesh.geometry.dim*mesh.geometry.dim,
                      x.shape[1]), dtype=np.float64)
    values[0] = x[0]
    values[1] = x[1]
    values[2] = x[0] + 2 * x[1]
    values[3] = x[1] - x[0]
    return values


form = dolfinx.fem.form(ufl.inner(F, F)*ufl.dx)
print(f"Pre {dolfinx.fem.assemble_scalar(form)}")
u.sub(1).interpolate(T)
print(f"Post {dolfinx.fem.assemble_scalar(form)}")

yielding

Pre 0.0
Post 3.5000000000000115
1 Like
5

Hey, thank you so much for your answers so far. What I did not completely understand is how you structured the function T(x):

  • Why do you square the mesh.geometry.dim?
  • Why do you assign only values[0] up to values[3], when values has a length of 9?
  • Which part of values exactly do I have to target to assign the identity matrix to F, i.e. the second function-part?

Thank you so much for your help, it really helps a lot!
Marco

6

As a Tensor element will have 3 * 3 entries if you use a cell with geometric/topological dimension 3.
The shape of values should be equal to the flattened version of the the tensor, i.e. values[0] = tensor[0,0], values[1] = tensor[0,1], values[2] = tensor[0,2], values[3] = tensor[1,0] etc.

Because this was just to illustrate to you how you should assign values for each component of the tensor.

Then you should create values with zeros, as I already showed

and set

values[0,:] = 1
values[3,:] = 1
values[6,:] = 1

You should verify this by writing F to file after interpolating data into it.

7

That seems to work for me. I cannot finally say until all other issues in my script are fixed, but according to printouts and such, it seems to do the job! Thank you very much for your time and help!