Dear FEniCS community,
I am working on a problem where my domain, i.e., the geometric model has to be parameterized. To tackle such a challenge, I am trying to map different versions/variants of this geometric model onto a reference geometry. Mapping function is a Laplace’s equation, so by solving a Laplace’s equation on reference domain, I get the new version of my domain. Now, the goal is to solve some fluid dynamics problem on this new version of domain, however, the idea is to formulate everything on reference domain to save the cost of meshing again and again for every new domain.
So, in order to formulate everything in reference domain, I need Jacobian matrix, which in this case is simply the gradient of solution of Laplace’s equation. Please have a look at the MWE to have an idea
from dolfin import *
geometry = Rectangle(Point(0., 0.), Point(4., 3.))
mesh = generate_mesh(geometry, 10)
class Top(SubDomain):
def inside(xself,x,on_boundary):
return x[1] > 3.0 - DOLFIN_EPS
class Bottom(SubDomain):
def inside(self,x,on_boundary):
return x[1] < DOLFIN_EPS
class Left(SubDomain):
def inside(xself,x,on_boundary):
return x[0] < DOLFIN_EPS
class Right(SubDomain):
def inside(xself,x,on_boundary):
return x[0] > 4.0 - DOLFIN_EPS
boundaries = MeshFunction("size_t", mesh, 1)
boundaries.set_all(0)
Left().mark(boundaries, 0)
Right().mark(boundaries, 1)
Top().mark(boundaries, 3)
Bottom().mark(boundaries, 2)
n = FacetNormal(mesh)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
u_trial = TrialFunction(V)
v_test = TestFunction(V)
exp = Expression(('0','0'),element=V.ufl_element())
g = exp
# Define variational problem
a = inner(grad(u_trial), grad(v_test))*dx
L = inner(g,v_test)*dx
# Boundary Conditions
bc = [(DirichletBC(V,Constant((0,0)),boundaries,0)),(DirichletBC(V,Constant((1,0)),boundaries,1))]
# Compute solution
sol = Function(V)
solve(a == L, sol, bc, solver_parameters={'linear_solver':'mumps'})
xderivative = sol.dx(0)
yderivative = sol.dx(1)
duoverdx = project(xderivative,V)
duoverdy = project(yderivative,V)
J = PETSc.Mat().createDense((len(u_trial),len(sol.vector())),array=([duoverdx.vector().vec().getArray()],[duoverdy.vector().vec().getArray()])).transpose()
This was my approach for obtaining the Jacobian matrix, however, the problem with this approach is that it is not a square matrix. I believe that the values in the matrix are correct, but this matrix is not sparse and square. What I could do is, I could manually make a square matrix of zeros and assign these values at the right positions according to the node numbers. But, this approach would very cumbersome and inefficient.
So, my question is, is there an efficient way to obtain the Jacobian matrix of such a system?
Thanks!
