Hello sir
I attached a screenshot, I want to know How we can compute G and g in FeniCS? The metric tensors correspond to the mappings to the reference element in FEM, and this is supported by UFL .
I think this form tis available in FEniCS. Please help me in this regard.
Consider the following:
import dolfinx
from mpi4py import MPI
import numpy as np
import ufl
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
K = ufl.JacobianInverse(mesh)
i, j, k = ufl.indices(3)
G = ufl.as_tensor(K[k, i]*K[k, j], (i, j))
g = ufl.as_vector(K[k, i], i)
Legacy DOLFIN
from dolfin import *
from ufl import JacobianInverse, indices
mesh = UnitSquareMesh(10, 10)
K = JacobianInverse(mesh)
i, j, k = indices(3)
G = as_tensor(K[k, i]*K[k, j], (i, j))
g = as_vector(K[k, i], i)
When I used above G and g in writing T_C and T_M which is given above screenshot then T_C shows an error.
T_M = pow(((1/(dt*dt))+inner(u,dot(u,G))+30*nu*nu*inner(G,G)),-0.5)
T_C = pow((dot(inner(T_M,g),g)),-1)
error
wloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices
(warning: it would break the library ABI, don't enable unless really needed).
hwloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices
(warning: it would break the library ABI, don't enable unless really needed).
Shapes do not match: <Power id=140670389092480> and <ComponentTensor id=140670388809424>.
Traceback (most recent call last):
File "/home/ayush/bb.py", line 36, in <module>
T_C = pow((dot(inner(T_M,g),g)),-1)
File "/usr/lib/python3/dist-packages/ufl/operators.py", line 158, in inner
return Inner(a, b)
File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 147, in __new__
error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: <Power id=140670389092480> and <ComponentTensor id=140670388809424>.
For me to help you you have to make a minimal reproducible example.
Minimal example
from dolfin import *
from ufl import JacobianInverse, indices
T = 1
num_steps = 10
dt = T/num_steps
mesh = UnitSquareMesh(2,2)
# Define Function Spaces
V = VectorElement("Lagrange", mesh.ufl_cell(), 1)
W = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
me = MixedElement([V, W, R])
Z= FunctionSpace(mesh,me)
###Define Test and trial function
up = Function(Z)
u, p, rho = split(up)
v, q, lamda = TestFunctions(Z)
u_exact = Expression(("20*x[0]*x[0]*(1-x[0])*(1-x[0])*x[1]*(2-3*x[1])","-20*x[1]*x[1]*(1-x[1])*(2*x[0]-6*x[0]*x[0]+4*x[0]*x[0]*x[0])"), domain=mesh, degree = 4, t=0)
u_n = project(u_exact,Z.sub(0).collapse())
# define weak formulation
nu=0.01
K = JacobianInverse(mesh)
i, j, k = indices(3)
G = as_tensor(K[k, i]*K[k, j], (i, j))
g = as_vector(K[k, i], i)
T_M = pow(((1/(dt*dt))+inner(u,dot(u,G))+30*nu*nu*inner(G,G)),-0.5)
T_C = pow((dot(inner(T_M,g),g)),-1)
error:
hwloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices
(warning: it would break the library ABI, don't enable unless really needed).
hwloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices
(warning: it would break the library ABI, don't enable unless really needed).
Shapes do not match: <Power id=140559681360192> and <ComponentTensor id=140559685521376>.
Traceback (most recent call last):
File "/home/ayush/cc.py", line 31, in <module>
T_C = pow((dot(inner(T_M,g),g)),-1)
File "/usr/lib/python3/dist-packages/ufl/operators.py", line 158, in inner
return Inner(a, b)
File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 147, in __new__
error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: <Power id=140559681360192> and <ComponentTensor id=140559685521376>.
I am trying to write equation 3.27 given in above screenshot.
T_M is a scalar value, so there is no way you can do an inner product with g,
I think you want
pow(T_M*dot(g, g), -1)
no it does not work
error:
hwloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices
(warning: it would break the library ABI, don't enable unless really needed).
hwloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices
(warning: it would break the library ABI, don't enable unless really needed).
Not expecting repeated indices.
Traceback (most recent call last):
File "/home/ayush/cc.py", line 31, in <module>
T_C = pow(T_M*dot(g,g),-1)
File "/usr/lib/python3/dist-packages/ufl/operators.py", line 176, in dot
return Dot(a, b)
File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 207, in __init__
fi, fid = merge_nonoverlapping_indices(a, b)
File "/usr/lib/python3/dist-packages/ufl/index_combination_utils.py", line 198, in merge_nonoverlapping_indices
error("Not expecting repeated indices.")
File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Not expecting repeated indices.
Some hard-coding of the loops should do it:
g = as_vector( [sum(K[k, i] for k in range(mesh.geometry().dim())) for i in range(mesh.geometry().dim())])
T_M = pow(((1/(dt*dt))+inner(u,dot(u,G))+30*nu*nu*inner(G,G)),-0.5)
T_C = pow((T_M*dot(g,g)),-1)
G = as_tensor(K[k, i]*K[k, j], (i, j))
g = as_vector(K[k, i], i)
Now, the code works. I have a doubt as you told in g ‘i’ and ‘k’ varies upto the dim of domain. No need to write such type of thing for G.
Thankyou so much
You need to write this for g as you do a negative power of it, and division does not support free indices in ufl. As G is in the nominator, it can have free indices (which is equivalent to summing over the dim of the domain).
sir, In T_M, we also take negative power (0.5) of G . I just ask this question for better understanding as i am little bit confused.
I’ve not used ufl indices for much in the past, so I cannot be more concrete as to describing what goes wrong.
If you do not like the usage of indices in G you can just hard-code it as I did for g.