Problem with creating vector with scalar function components

I will go over the suggested links. However, in the meantime, I am facing another problem with another script. In this one, c (x, y) is a scalar function and the domain is an annulus where the inner circle is replaced by a random closed curve. The inner curve has the physical tag 1 and the outer circle has the physical tag 2. Then the variational form is

W = FunctionSpace(mesh, "Lagrange", 1)
w = TestFunction(W)
c = TrialFunction(W)
c = Function(W)
a = dot(grad(c), grad(w))*dx
L = -((((k/D)*c)*w)*dx) + inner(grad(c)*n, w)*ds_custom(1) + inner(grad(c)*n, w)*ds_custom(2)

This is returning the error

Invalid ranks 1 and 1 in product.

Since grad(c) is a vector I tried replacing W = VectorFunctionSpace(mesh, "Lagrange", 1) and that gives the following error:

Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2, 2) and free indices with labels ().

What am I doing wrong here?

First of all, you need to supply a minimal reproducible code that one can run after copy-paste and being able to reproduce the error. Please also think carefully about what operations you are trying to apply.
Following is a minimal code example I have interpreted from your code above, that does not produce an error. Note that you are overwriting c in your code, and therefore I have changed the formulation to a residual formulation using a Function, such that you can solve the problem with solve(F==0, c, bcs=...).

Please read carefully through the following code and make sure you understand it.

from dolfin import *
mesh = UnitSquareMesh(10, 10)
W = FunctionSpace(mesh, "Lagrange", 1)
k = Constant(0.1)
D = Constant(0.2)
w = TestFunction(W)
c = TrialFunction(W)
c = Function(W)
n = FacetNormal(mesh)

F = inner(grad(c), grad(w))*dx
F -= -(k/D)*c*w*dx + dot(grad(c), n)*w*ds + dot(grad(c), n)*w*ds
1 Like

My code worked but I didn’t quite get the last line

F -= -(k/D)*c*w*dx + dot(grad(c), n)*w*ds + dot(grad(c), n)*w*ds

Does the last line mean

F = inner(grad(c), grad(w))*dx - (k/D)*c*w*dx + dot(grad(c), n)*w*ds + dot(grad(c), n)*w*ds ?

Your original equation is: a==L where

I have now rewritten this as F=a-L=0
thus the sign of every term in the last line is opposite of what you have stated, as it is subtracting each term.

I got that! Thanks for the explanation. But I have a question. Why doesn’t the following code work? Is it because we are dealing with a function space and not a vector function space? Whenever I am trying to execute the below code, it’s returning the error message

    Error:   Unable to define linear variational problem a(u, v) == L(v) for all v.
    *** Reason:  Expecting the left-hand side to be a bilinear form (not rank 1).

Code:

a = inner(grad(c), grad(w))*dx
L =  -((((k/D)*c)*w)*dx) - inner(gamma1, w)*ds_custom(1)

# Compute solution
solve(a == L, c, bcs)

You need to define c as a TrialFunction if you want to solve something on the form a==L

But I did define c as c = TrialFunction(W). Below is the entire code along with the gmsh script of the geometry.

Python script:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt


from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp



# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

import meshio
msh = meshio.read("intimapdgf.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)

# Normal vector to the intima
n = FacetNormal(mesh)

# Create mesh and define function space

W = FunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition and neumann boundary functions
D, k, t = 0.7, 0.3, 3
f0, ts = 5, 2
gamma1 = (-1/D)*Constant(f0*t/(ts + t)) # Neumann condition on inner curve -dot(grad(c),n)=gamma1
gamma2 = Constant(0.0)    # Dirichlet condition on outer circle 

bc2 = DirichletBC(W, gamma2, mf, 2)

bcs=[bc2]

# Define variational problem
w = TestFunction(W)
c = TrialFunction(W)
c = Function(W)
a = inner(grad(c), grad(w))*dx
a -= -((((k/D)*c)*w)*dx) - inner(gamma1, w)*ds_custom(1)


# Compute solution
solve(a == 0, c, bcs)


# Save solution to file in PVD format for Paraview
file = File("intimapdgf.pvd");
file << c;

gmsh script:

//+
SetFactory("OpenCASCADE");
//+
Circle(2) = {0, 0, 0, 3, 0, 2*Pi};
//+
Point(3) = {-1.6, 0.3, 0, 1.0};
//+
Point(4) = {-1.4, -0.7, 0, 1.0};
//+
Point(5) = {-1.1, -1.3, 0, 1.0};
//+
Point(6) = {-0.3, -1.7, 0, 1.0};
//+
Point(7) = {0.8, -2.1, 0, 1.0};
//+
Point(8) = {1.7, -1.3, 0, 1.0};
//+
Point(9) = {1.7, -0.4, 0, 1.0};
//+
Point(10) = {1.9, 0.5, 0, 1.0};
//+
Point(11) = {1.6, 1.2, 0, 1.0};
//+
Point(12) = {0.9, 1.5, 0, 1.0};
//+
Point(13) = {0, 1.5, 0, 1.0};
//+
Point(14) = {-0.6, 0.9, 0, 1.0};
//+
Spline(3) = {3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 3};
//+
Physical Curve("C2", 2) = {2};
//+
Physical Curve("C3", 3) = {1};
//+
Hide "*";
//+
Show {
  Point{3}; Curve{3}; 
}
//+
Curve Loop(1) = {2};
//+
Curve Loop(2) = {3};
//+
Plane Surface(1) = {1, 2};
//+
Hide "*";
//+
Show {
  Point{1}; Point{2}; Point{3}; Curve{1}; Curve{2}; Curve{3}; 
}
//+
Curve Loop(3) = {1};
//+
Curve Loop(4) = {2};
//+
Plane Surface(2) = {3, 4};
//+
Show "*";
//+
Physical Surface("intima", 4) = {1};
//+
Hide "*";
//+
Show {
  Point{1}; Point{2}; Point{3}; Curve{1}; Curve{2}; Curve{3}; Surface{1}; Surface{2}; 
}
//+
Physical Curve("C1", 1) = {3};

You are overwriting c in the following line. You Need to remove the second line

I tried that. That gives the following error with the below modification of the script

w = TestFunction(W)
c = TrialFunction(W)
a = inner(grad(c), grad(w))*dx
L = -((((k/D)*c)*w)*dx) - inner(gamma1, w)*ds_custom(1)


# Compute solution
solve(a == L, c, bcs)

Error:

Expecting second argument to be a Function.

You should approach this as:

w = TestFunction(W)
c = TrialFunction(W)
a = inner(grad(c), grad(w))*dx
L = -((((k/D)*c)*w)*dx) - inner(gamma1, w)*ds_custom(1)

ch = Function(W)
# Compute solution
solve(a == L, ch, bcs)

The solution of your problem is then contained in the variable ch

Note that in my post above i didn’t realize that you have moved one of the terms with c back to the RHs. This should be moved to the left hand side (a)

Ok. But why do I need to introduce a new variable ch? Why can’t the solution be assigned to c itself? Below is the part of a code which worked perfectly fine.

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

Here I am overwriting u. So why the above works perfectly fine whereas the below doesn’t?

w = TestFunction(W)
c = TrialFunction(W)
c = Function(W)
a = inner(grad(c), grad(w))*dx
L = -((((k/D)*c)*w)*dx) - inner(gamma1, w)*ds_custom(1)


# Compute solution
solve(a == L, c, bcs)

You can use c for both cases, But it is all about the code order. You have to define the variational form a and L with the trial function before redefining c as a function. In general in think it is bad coding pratice to redefine variables, and therefore I use ch and not c to define the variable containing the numerical solution.

Ok, got that! Thank you so much for the explanation.

I want to discuss another problem with you. It is related to all the problems we discussed in this post so far. I have solved three separate problems, one for u, one for Sn[2] and the last one for ch, which we just discussed. All the problems are solved using different python scripts on different subdomains of a composite domain. When I am trying to combine all of them, it’s giving me an empty solution. I can’t figure out the problem. I read some post regarding solving variational problems on different subdomains in the same script. But I didn’t understand most of those. It would be great if I can discuss it with you. I am still trying to figure it out by myself. If I feel totally stuck, I will get back to you.

Thank you for your time and patience.

You need to sit down and carefully think about how to formulate your problem as a minimal reproducable example.

1 Like

Hi dokken! Here is a toy problem I created. Below are the gmsh script (named as example.msh) followed by the python script. In the problem, the composite domain is a disk and an annulus. Problem 1 is defined on the disk and problem 2 is defined on the annulus (I have marked that in the script). I used MeshView.create to extract the subdomains, defined the measure on the domain. But I am unable to define a Dirichlet boundary condition on curve 1.

gmsh script:

//+
SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 0, 0, 5, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 3, 0, 2*Pi};
//+
Curve Loop(1) = {2};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {1};
//+
Curve Loop(3) = {2};
//+
Plane Surface(2) = {2, 3};
//+
Physical Curve(1) = {2};
//+
Physical Curve(2) = {1};
//+
Physical Surface(5) = {1};
//+
Physical Surface(4) = {2};

Python script:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt


from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp



# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

import meshio
msh = meshio.read("example.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)


s5 = MeshView.create(mf, 5)
s4 = MeshView.create(mf, 4)

W1 = FunctionSpace(s5, "Lagrange", 1)
W2 = FunctionSpace(s4, "Lagrange", 1)

V = MixedFunctionSpace(W1, W2)

dx5 = Measure("dx", domain=s5)
dx4 = Measure("dx", domain=s4)

# PROBLEM FOR SURFACE 5

# Normal vector to s5
n5 = FacetNormal(s5)

# Create mesh and define function space

# Define boundary condition
G , mu = 1, 0.1
u_D=Constant(0.0)

bc = DirichletBC(V.sub_space(0), u_D, mf, 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx5
L = -f*v*dx5

# Compute solution
us = Function(V)
solve(a == L, us, bc)

# PROBLEM FOR SURFACE 4

D, k, t = 0.7, 0.3, 3
f0, f1, beta = 5, 2, 0.3
gamma1 = (-1/D)*Constant(f0 + f1*(1-exp(-beta*(t**2)))) # Neumann condition on inner curve -dot(grad(c),n)=gamma1
gamma2 = Constant(0.0)    # Dirichlet condition on outer circle 

bc2 = DirichletBC(V.sub_space(1), gamma2, mf, 2)

bcs=[bc2]

# Define variational problem
w = TestFunction(V)
c = TrialFunction(V)
a = inner(grad(c), grad(w))*dx4 + (((k/D)*c)*w)*dx4
a -= - inner(gamma1, w)*ds_custom(1)
 

# Compute solution
cs= Function(V)
solve(a == 0, cs, bcs)


# Save solution to file in PVD format for Paraview
file = File("example.pvd");
file << cs;

The error is as follows:

Error:   Unable to create Dirichlet boundary condition.
*** Reason:  User MeshFunction and FunctionSpace meshes are different.

Where am I wrong?

Hi @dokken , following the example demo_meshview-2D2D.py, I have modified the above python script as follows (the geometry is the same as in the gmsh script posted before):

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt


from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp



# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

import meshio
msh = meshio.read("example.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)


s5 = MeshView.create(cf, 5) # Extracting subspace 5
s4 = MeshView.create(cf, 4) # Extracting subspace 4

W1 = FunctionSpace(s5,"Lagrange",1) # Defining function space on s5
W2 = FunctionSpace(s4,"Lagrange",1) # Defining function space on s4

def bs5(x):  # Defining boundary for subspace s5, curve 1
  return mf == 1 
  
def b2s4(x): # Defining boundary for subspace s4, curve 2
  return mf == 2
  
# PROBLEM FOR DOMAIN 5 WITH BOUNDARY 1 
  
# Define boundary condition
G , mu = 1, 0.1
u_D=Constant(0.0)

bc = DirichletBC(W1, u_D, bs5)

# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx

# Compute solution
us = Function(W1)
solve(a == L, us, bc)

# Save solution to file in PVD format for Paraview
file = File("examplelumenvelocity.pvd");
file << us;

# PROBLEM FOR SURFACE 4 WITH BOUNDARIES 1 & 2

D, k, t = 0.7, 0.3, 3
f0, f1, beta = 5, 2, 0.3
gamma1 = (-1/D)*Constant(f0 + f1*(1-exp(-beta*(t**2)))) # Neumann condition on inner curve -dot(grad(c),n)=gamma1
gamma2 = Constant(0.0)    # Dirichlet condition on outer circle 

bc2 = DirichletBC(W2, gamma2, b2s4)

bcs=[bc2]

# Define variational problem
w = TestFunction(W2)
c = TrialFunction(W2)
a = inner(grad(c), grad(w))*dx + (((k/D)*c)*w)*dx
a -= - inner(gamma1, w)*ds_custom(bs5)
 

# Compute solution
cs= Function(W2)
solve(a == 0, cs, bcs)


# Save solution to file in PVD format for Paraview
file = File("example.pvd");
file << cs;

This returns the error:

*** Warning: Found no facets matching domain for boundary condition.
Invalid subdomain_id <function bs5 at 0x7fc435e3f1f0>.

How to fix this?

Hi @dokken, there seems to be a problem in the way I am defining the boundaries to each of the submeshes, i.e.,

def bs5(x):  # Defining boundary for subspace s5, curve 1
  return mf == 1 
  
def b2s4(x): # Defining boundary for subspace s4, curve 2
  return mf == 2

where 1 and 2 are the physical tags attached to the boundaries of the parent mesh. What am I missing here?

Please stop adding to this thread with new questions. Your question is now on a vastly different topic than what this thread started with. Please start a new thread and reduce your problem to a minimal reproducible code using a built-in mesh to illustrate the issues.

1 Like

Hi @dokken , my apologies for posting on the same thread. I am sorry for that.

I had already created a new thread yesterday [(https://fenicsproject.discourse.group/t/how-to-define-bcs-on-boundaries-of-submeshes-of-a-parent-mesh/5470)], but haven’t received any response from anyone yet. So I thought of continuing this thread.

Since my problem concerns with extraction of physical tags of boundaries from gmsh to FEniCS when defining subdomains and defining BCs on them, I am not sure if a built-in mesh would be an appropriate example to address the problem. In the link mentioned above, I have tried my best to reduce the problem to a minimal working example to illustrate the issue with a mesh imported from gmsh.

Again, I am sorry for continuing this thread. My sincere apologies for any inconvenience this may have caused.