# Discrepancy, with assemble and assmbled operator

**URL:** <https://fenicsproject.discourse.group/t/discrepancy-with-assemble-and-assmbled-operator/5057>\
**Category:** mathematics\
**Created:** [February 9, 2021, 1:25am UTC](https://fenicsproject.discourse.group/t/discrepancy-with-assemble-and-assmbled-operator/5057 "2021-02-09T01:25:19Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![tucker](https://avatars.discourse-cdn.com/v4/letter/t/8c91f0/32.png) [@tucker](https://fenicsproject.discourse.group/u/tucker)\
**Post date:** [February 9, 2021, 1:25am UTC](https://fenicsproject.discourse.group/t/discrepancy-with-assemble-and-assmbled-operator/5057/1 "2021-02-09T01:25:19Z")

</div>

I have created a mesh in `gmsh` for the Haut Glacier d’Arolla ([TC - Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP–HOM)](http://www.the-cryosphere.net/2/95/2008/)) and converted to `.xml` format for use in `FEniCS` using `dolfin-convert`.

I decided to test my mesh by solving the following Poisson problem.

 ![Screenshot from 2021-02-08 16-51-31](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/de0e660806d4e0860bd9e8bb592b8b1345fa4a7d.png)

After solving the discretized problem, checking that the solution is reasonable by visualizing in `ParaView` I then check the H1 semi-norm of the solution, defined as

 ![Screenshot from 2021-02-08 17-07-24](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/8/80276afb2979930596e94cf57ef78a4268e08df3.png)

in two ways. First by using the `assemble` command in `dolfin`, the numerical value returned is `nan`. There is an alternative way to compute the H1 semi-norm. That is

 ![Screenshot from 2021-02-08 17-14-16](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/8/8ef4d50d50ad38ddd727224b90c673a4522a17f2.png)

Where **A** is the finite element stiffness matrix. **u** ^T **A u** evaluates to a reasonable nonnegative finite number.

I have included the full code below., but of course you will not be able to test without the mesh.

I am hoping to get insight into why this would be the case, what goes on internally in `FEniCS` when one assembles a finite element system matrix, such as **A** , with Dirichlet conditions that causes the discrepancy that I am seeing.

> import dolfin as dl  
> import numpy as np

```
# mesh of PDE domain
meshdir = "../"
meshname = "Arolla2"
mesh = dl.Mesh(meshdir+meshname+".xml")
boundary_markers = dl.MeshFunction("size_t", mesh,\
                   meshdir+meshname+"_facet_region.xml")

# discrete function spaces
P1 = dl.FiniteElement("CG", mesh.ufl_cell(), 1)
Vh = dl.FunctionSpace(mesh, P1)

# Homogeneous Dirichlet conditions
bc = [dl.DirichletBC(Vh, dl.Constant(0.0), boundary_markers, 1),\
      dl.DirichletBC(Vh, dl.Constant(0.0), boundary_markers, 2)]

# Right hand side forcing term
f = dl.Constant(1.0)

vh = dl.TestFunction(Vh)
uh = dl.TrialFunction(Vh)
a = dl.inner(dl.grad(uh), dl.grad(vh))*dl.dx
L = f*vh*dl.dx
A, b = dl.assemble_system(a, L, bc)

Asolver = dl.PETScKrylovSolver()
dl.PETScOptions.set("ksp_monitor")
dl.PETScOptions.set("ksp_monitor_true_residual")
dl.PETScOptions.set("ksp_type", "cg")
dl.PETScOptions.set("ksp_rtol", "1.e-12")
dl.PETScOptions.set("ksp_atol", "1.e-12")
Asolver.set_from_options()
Asolver.set_operator(A)

u = dl.Function(Vh)
Asolver.solve(u.vector(), b)

# for the determined solution,
# compute sqrt(int_(Omega) ||grad(u)||^2 dV) 
u_H1norm_1 = np.sqrt(dl.assemble(dl.inner(dl.grad(u), dl.grad(u))*dl.dx(mesh))) 
print("||u||_H^1(Omega)^2 = %1.3e (directly assembled)" %(u_H1norm_1))

# vector to hold the product of the solution with
# the finite element system matrix A
Au = dl.Vector()
A.init_vector(Au, 0)

# Au = A*u
A.mult(u.vector(), Au)

# sqrt(u^T Au)
u_H1norm_2 = np.sqrt(Au.inner(u.vector()))
print("||u||_H^1(Omega) = %1.3e (from assembled operators)" \
             %(u_H1norm_2))

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 9, 2021, 8:56am UTC](https://fenicsproject.discourse.group/t/discrepancy-with-assemble-and-assmbled-operator/5057/2 "2021-02-09T08:56:08Z")

</div>

I cannot reproduce this with a built in mesh:

```auto
import dolfin as dl
import numpy as np

mesh = dl.UnitSquareMesh(30,30)
# discrete function spaces
P1 = dl.FiniteElement("CG", mesh.ufl_cell(), 1)
Vh = dl.FunctionSpace(mesh, P1)

# Homogeneous Dirichlet conditions
bc = [dl.DirichletBC(Vh, dl.Constant(0.0), "on_boundary")]

# Right hand side forcing term
f = dl.Constant(1.0)

vh = dl.TestFunction(Vh)
uh = dl.TrialFunction(Vh)
a = dl.inner(dl.grad(uh), dl.grad(vh))*dl.dx
L = f*vh*dl.dx
A, b = dl.assemble_system(a, L, bc)

Asolver = dl.PETScKrylovSolver()
dl.PETScOptions.set("ksp_monitor")
dl.PETScOptions.set("ksp_monitor_true_residual")
dl.PETScOptions.set("ksp_type", "cg")
dl.PETScOptions.set("ksp_rtol", "1.e-12")
dl.PETScOptions.set("ksp_atol", "1.e-12")
Asolver.set_from_options()
Asolver.set_operator(A)

u = dl.Function(Vh)
Asolver.solve(u.vector(), b)

# for the determined solution,
# compute sqrt(int_(Omega) ||grad(u)||^2 dV) 
u_H1norm_1 = np.sqrt(dl.assemble(dl.inner(dl.grad(u), dl.grad(u))*dl.dx(mesh))) 
print("||u||_H^1(Omega)^2 = %1.3e (directly assembled)" %(u_H1norm_1))

# vector to hold the product of the solution with
# the finite element system matrix A
Au = dl.Vector()
A.init_vector(Au, 0)

# Au = A*u
A.mult(u.vector(), Au)

# sqrt(u^T Au)
u_H1norm_2 = np.sqrt(Au.inner(u.vector()))
print("||u||_H^1(Omega) = %1.3e (from assembled operators)" \
             %(u_H1norm_2))

```

which returns:

```auto
||u||_H^1(Omega)^2 = 1.871e-01 (directly assembled)
||u||_H^1(Omega) = 1.871e-01 (from assembled operators)

```

(using docker and: (`docker run -it -v $(pwd):/home/fenics/shared -w /home/fenics/shared --rm quay.io/fenicsproject/dev:latest`)  
Therefore, it is most likely your mesh that has issues.  
I also reproduced it with an unstructured mesh, using mshr:

```auto
cylinder = mshr.Cylinder(dl.Point(0, 0, 0), dl.Point(0, 0, -9), 8.0, 8.0)
domain = cylinder
mesh = mshr.generate_mesh(domain, 20)

```

with result:

```auto
||u||_H^1(Omega)^2 = 7.148e+01 (directly assembled)
||u||_H^1(Omega) = 7.148e+01 (from assembled operators)

```

---

<div class="post-metadata">

**Author:** ![tucker](https://avatars.discourse-cdn.com/v4/letter/t/8c91f0/32.png) [@tucker](https://fenicsproject.discourse.group/u/tucker)\
**Post date:** [February 12, 2021, 5:59am UTC](https://fenicsproject.discourse.group/t/discrepancy-with-assemble-and-assmbled-operator/5057/3 "2021-02-12T05:59:32Z")

</div>

The issue appears to have been due to an issue with how I generated the mesh in gmsh. I was still hoping to get insight into the specifics of how a matrix is altered when Dirichlet conditions are applied, as I expected the results above to be the same.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 12, 2021, 6:43am UTC](https://fenicsproject.discourse.group/t/discrepancy-with-assemble-and-assmbled-operator/5057/4 "2021-02-12T06:43:15Z")

</div>

To assemble a matrix system symmetrically with Dirichlet conditions, the most common thing to do is to apply lifting of the boundary conditions. This is done by ignoring all contributions from cells containing Dirichlet-dofs during assembly, and subtract the matrix contribution for the Dirichlet dofs from the RHS.  
Each row or column corresponding to a Dirichlet condition will simply have one on the diagonal, and be zero everywhere else.

---

<div class="post-metadata">

**Author:** ![tucker](https://avatars.discourse-cdn.com/v4/letter/t/8c91f0/32.png) [@tucker](https://fenicsproject.discourse.group/u/tucker)\
**Post date:** [February 12, 2021, 5:08pm UTC](https://fenicsproject.discourse.group/t/discrepancy-with-assemble-and-assmbled-operator/5057/5 "2021-02-12T17:08:35Z")

</div>

Thank you so much, this is helpful.