Real and Imaginary part of solution vector - dolfinx

1

Hi,
I am currently performing a harmonic analysis of the damped 3D elastic model in FenicsX. I have implemented the following code below. I am stuck in splitting the solution vector into its real and imaginary parts. I want to plot the frequency response v/s amplitude values of the z-component of the displacement field. In order to do so, I want to access the real and imaginary values of the displacement fields z-component.
Amplitude = sqrt(Uz_real^^2 +Uz_imag^^2).
In this case, how can I access the real and imaginary values in dolfinx? rather than exporting the result to xdmf file to view in ParaView. Is there a way to split the solution into real and imaginary parts?

import dolfinx
from dolfinx.fem import (Constant,dirichletbc, Function, FunctionSpace, VectorFunctionSpace,LinearProblem,assemble_scalar, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem import (assemble_matrix, assemble,assemble_vector, form, apply_lifting, locate_dofs_geometrical, set_bc)
import ufl
from ufl import (Form, SpatialCoordinate, TestFunction, TrialFunction,dx, ds, grad, FacetNormal,inner, max_value,nabla_grad, nabla_div, Identity)
from mpi4py import MPI
import dolfinx.io
import numpy as np
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from dolfinx.io import XDMFFile


with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
    m_subdomains = xdmf.read_meshtags(mesh, name="Grid")

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mf.xdmf", "r") as xdmf:
    m_boundaries = xdmf.read_meshtags(mesh, name="Grid")

# function space
V = VectorFunctionSpace(mesh, ("CG", 2))
num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
print(f"Number of dofs global: {num_dofs_global}")
dx = ufl.Measure("dx", domain=mesh, subdomain_data=m_subdomains)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=m_boundaries)

# Material parameters and elastic constitutive relations 
E, nu = (2E11), (0.31)
rho = (7900)

# Lame coefficient for constitutive relation
mu = E/2./(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

# Strain function
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

# Stress function
def sigma(u):
    return lmbda * nabla_div(u) * Identity(u.geometric_dimension()) + 2*mu*epsilon(u)

# Rayliegh Coefficients
alpha = (0) # mass coefficient
beta = (1e-3) # stiffness coefficient
# trial and test functions
u = TrialFunction(V)  
w = TestFunction(V)

# dirichlet bc
left_facets = m_boundaries.indices[m_boundaries.values==1]
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)

u_D = Function(V)
u_D = np.array([0,0,1e-3], dtype=ScalarType)
bc = [dirichletbc(u_D, left_dofs, V)]

T = Constant(mesh, ScalarType((0, 0, 0)))
start = 100; d_freq = 10;  freq=130
with XDMFFile(MPI.COMM_WORLD,"Frequency_Response.xdmf", "w") as outfile:
    outfile.write_mesh(mesh)
    while start <= freq:
        print(start)
        omega=2*np.pi*start
        print(omega)
        # Matrices
        M_form = rho*inner(u, w)*dx
        K_form = inner(sigma(u), epsilon(w))*dx
        C_form = alpha * M_form + beta * K_form
        a =  - omega**2*M_form  + 1j * omega * C_form + K_form
        L =  inner(T, w) * ds
        # Set up the PDE
        solnU = Function(V)
        problem = LinearProblem(a, L, bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
        solnU = problem.solve()
        outfile.write_function(solnU, freq)
        start = start + d_freq
    else:
        print("Harmonic Analysis Completed")
2

Consider the following MWE:

import dolfinx.mesh
import dolfinx.fem
from mpi4py import MPI
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 3, 3)
V = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))

u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0]+1j*x[1])


uR = dolfinx.fem.Function(V, dtype=np.float64)
uI = dolfinx.fem.Function(V, dtype=np.float64)

uR.x.array[:] = u.x.array.real
uI.x.array[:] = u.x.array.imag
print(u.x.array, "\n", uR.x.array, "\n", uI.x.array)

uAmp = dolfinx.fem.Function(V, dtype=np.float64)
uAmp.x.array[:] = np.sqrt(uR.x.array**2+uI.x.array**2)
print(uAmp.x.array)
2 Likes