problem with definning mixed product space


149
views
1
3 months ago by
Ali  
Hi All,

I am trying to print out stiffness matrix regards to the bilinear form below:
\begin{equation}
a( (\boldsymbol{Z} , \boldsymbol{Q}) , (\boldsymbol{\Upsilon} ,\boldsymbol{R} ) )= ((\boldsymbol{Q}_{1},\mathbf{grad}\,\boldsymbol{\Upsilon}_{1})) + ((\boldsymbol{Q}_{2},\mathbf{grad}\,\boldsymbol{\Upsilon}_{2})) + \\ ((\boldsymbol{Z}_{1},\boldsymbol{\Upsilon}_{1})) + (( \mathbf{grad}\boldsymbol{Z}_{1},\mathbf{grad}\,\boldsymbol{\Upsilon}_{1})) + ((\boldsymbol{Z}_{2},\boldsymbol{\Upsilon}_{2})) + (( \mathbf{grad}\boldsymbol{Z}_{2},\mathbf{grad}\,\boldsymbol{\Upsilon}_{2})) + \\
((\boldsymbol{Q}_{1},\boldsymbol{R}_{1})) + (( \mathbf{div}\boldsymbol{Q}_{1},\mathbf{div}\,\boldsymbol{R}_{1})) +
((\boldsymbol{Q}_{2},\boldsymbol{R}_{2})) + (( \mathbf{div}\boldsymbol{Q}_{2},\mathbf{div}\,\boldsymbol{R}_{2})) \\
( \boldsymbol{Z} ,  \boldsymbol{\Upsilon} ) \in H^1 \space and \space
( \boldsymbol{Q} ,  \boldsymbol{R}) \in H^d
\end{equation}

""((  ,  )) denotes both the L2 inner product of vector fields and tensor fields.""


Please find full code in below:

from dolfin import *
import scipy
from scipy.linalg import sqrtm
from petsc4py import PETSc
from numpy import array2string
import numpy as np
import sys
import time


start_time = time.time()
# Create mesh
mesh = UnitSquareMesh(2, 2)
# Define function spaces
RT = FiniteElement("RT", mesh.ufl_cell(), 1)
CG = FiniteElement("CG", mesh.ufl_cell(), 1)
W = (CG*CG) * (RT*RT)
V = FunctionSpace(mesh,W)
w = Function(V)
#Define DOFs
Gdof = V.dim()
# Define trial and test functions
(Z1 ,Q1) = TrialFunctions(V)
(Z2 ,Q2) = TrialFunctions(V)
(yup1 ,R1) = TestFunctions(V)
(yup2 ,R2) = TestFunctions(V)
####################################################
# Define bilinear form and right hand side

a = (dot(Q1, grad(yup1)) + dot(Q2, grad(yup2))  \
+ dot(Z1, yup1) + dot(Z2, yup2) + dot(grad(Z1), grad(yup1)) + dot(grad(Z2), grad(yup2)) \
+ dot(Q1, R1) + dot(Q2, R2) + ((div(Q1) * div(R1))) + ((div(Q2) * div(R2))))*dx

#####################################################

print('Global_DOF',Gdof)

matA = assemble(a)

viewer = PETSc.Viewer().createASCII("total01.m","w")
viewer.pushFormat(PETSc.Viewer.Format.ASCII_MATLAB)
viewer(as_backend_type(matA).mat())

​

The error is :
ufl.log.UFLException: Dimension mismatch in dot product.​

I have defined mixed finite element space as follows:
W = (CG*CG) * (RT*RT)​

That I think is a source of error but I do not know how can I fix that.


Thanks,
Ali


Community: FEniCS Project

1 Answer


2
3 months ago by
I am sorry, I don't really understand what you want to do exactly. But
  • inner must be used for the inner product between tensors
  • you cannot use two Test/Trial functions in the same form, if I understand it  $Z$Z  and $\Upsilon$ϒ  are vectors, the components of which can be accessed by Z[0] Z[1] for representing your  $Z_1$Z1 and  $Z_2$Z2  variables. In that case, you would need to define W as the product of VectorFunctionSpace
Thanks for your comment.

The goal of this script is finding stiffness matrix of the bilinear form a.

I have changed my code based on what you said but it is not working.

from dolfin import *
import scipy
from scipy.linalg import sqrtm
from petsc4py import PETSc
from numpy import array2string
import numpy as np
import sys
import time


start_time = time.time()
# Create mesh
mesh = UnitSquareMesh(2, 2)
# Define function spaces
RT = FiniteElement("RT", mesh.ufl_cell(), 1)
CG = FiniteElement("CG", mesh.ufl_cell(), 1)
W = CG*RT
V = VectorFunctionSpace(mesh,W)
w = Function(V)
#Define DOFs
Gdof = V.dim()
# Define trial and test functions
(Z ,Q) = TrialFunctions(V)
(yup ,R) = TestFunctions(V)

# Define bilinear form and right hand side

a = (dot(Q[0], grad(yup[0])) + dot(Q[1], grad(yup[1])) \
+ dot(Z[0], yup[0]) + dot(Z[1], yup[1]) + dot(grad(Z[0]), grad(yup[0])) + dot(grad(Z[1]), grad(yup[1])) \
+ dot(Q[0], R[0]) + dot(Q[1], R[1]) + ((div(Q[0]) * div(R[0]))) + ((div(Q[1]) * div(R[1]))))*dx

#####################################################

print('Global_DOF',Gdof)

matA = assemble(a)

viewer = PETSc.Viewer().createASCII("total01.m","w")
viewer.pushFormat(PETSc.Viewer.Format.ASCII_MATLAB)
viewer(as_backend_type(matA).mat())​

The Error is:

TypeError: VectorFunctionSpace() takes at least 3 arguments (2 given)
written 3 months ago by Ali  
That is not what I was suggesting. Change to
RT = VectorElement("RT", mesh.ufl_cell(), 1, dim=2)
CG = VectorElement("CG", mesh.ufl_cell(), 1, dim=2)
W = CG*RT
V = FunctionSpace(mesh,W)​
written 3 months ago by bleyerj  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »