# Unable to successfully call PETSc function 'MatSetValuesLocal'

**URL:** <https://fenicsproject.discourse.group/t/unable-to-successfully-call-petsc-function-matsetvalueslocal/8172>\
**Category:** variational formulation\
**Created:** [April 27, 2022, 8:40pm UTC](https://fenicsproject.discourse.group/t/unable-to-successfully-call-petsc-function-matsetvalueslocal/8172 "2022-04-27T20:40:51Z")\
**Posts on this page:** 1\
**Showing post:** 5

<div class="post-metadata">

**Author:** ![bbataszai](https://avatars.discourse-cdn.com/v4/letter/b/8c91f0/32.png) [@bbataszai](https://fenicsproject.discourse.group/u/bbataszai)\
**Post date:** [May 3, 2022, 9:10am UTC](https://fenicsproject.discourse.group/t/unable-to-successfully-call-petsc-function-matsetvalueslocal/8172/5 "2022-05-03T09:10:28Z")

</div>

Hi Cécile,

Thank you a lot. Oh my god, I checked the code so often to not have a stupid error like this in it. Sorry for bothering you with this.

Something I don’t understand happens with the measures: Depending on how they are specified to integrate the last term of `a2` only over the domain T the code either results in a ` Error: Unable to successfully call PETSc function 'MatSetValuesLocal'` or it solves. Can you help me understand why this is the case? (perhaps related to [this](https://fenicsproject.discourse.group/t/islocaltoglobalmappingapply-error/1260/2) post from 2019)

Full, running examples are below. To result in the error I specify restrict the \int\_T \text{d}x M\_t \phi\_{et} \psi\_{HaT} term in `a2` as follows:

```python
# define dx_torso with subdomain data
dx_heart = df.Measure("dx", mesh_heart)
dx_torso = df.Measure("dx", mesh, subdomain_data=subdomains)
# Code as before [...]
a2 = (
  M_i * df.inner(df.grad(v_mt), df.grad(psi_HaT)) * dx_heart 
  + (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_heart
)
# Integrate only over $T$ -- next line: dx_torso -> dx_torso(0)
a2+= M_t * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_torso(0)

```

To obtain running code, I use the measure that goes with the trial function:

```python
dx_heart = df.Measure("dx", mesh_heart)
dx_torso = df.Measure("dx", mesh, subdomain_data=subdomains)
# Code as before [...]

# Measures in line 
a2 = (
   M_i * df.inner(df.grad(v_mt), df.grad(psi_HaT)) * dx_heart
   + (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HaT)) *dx_torso(1)
   + M_t * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_torso(0)
)

```

It is kind of difficult to ask a clear question… Does this mean I have to be extra careful, which measure I use? Is the rule of thumb, to use the measure that goes with the trial function a good one?

The example codes are here:

> **MatSetValuesLocal - error example:**
>
> ```python
> import dolfin as df
> 
> mesh = df.UnitSquareMesh(20, 20)
> interior = df.CompiledSubDomain("std::max(fabs(x[0] - 0.5), fabs(x[1] - 0.5)) < 0.25")
> subdomains = df.MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
> # mark domain
> for cell in df.cells(mesh):
> x = cell.midpoint().array()
> subdomains[cell] = int(interior.inside(x, False))
> 
> # create interior mesh
> mesh_heart = df.MeshView.create(subdomains, 1)
> 
> # Function spaces
> H = df.FunctionSpace(mesh_heart, "Lagrange", 1)
> HaT = df.FunctionSpace(mesh, "Lagrange", 1)
> W = df.MixedFunctionSpace(H, HaT)
> 
> # Variational form
> dx_heart = df.Measure("dx", mesh_heart)
> dx_all = df.Measure("dx", mesh, subdomain_data=subdomains)
> 
> psi_H, psi_HaT = df.TestFunctions(W)
> v_mt, phi_et = df.TrialFunctions(W)
> 
> # last timestep v and ionic current
> # can this be done like this?
> v_mtm1 = df.Function(H)
> I_ion = df.Function(H)
> 
> # conductivities
> M_i = df.Constant(1)
> M_e = df.Constant(1)
> M_t = df.Constant(1)
> 
> # eq1
> a1 = (
> M_i * df.inner(df.grad(v_mt), df.grad(psi_H)) * dx_heart
> + M_i * df.inner(df.grad(phi_et), df.grad(psi_H)) * dx_heart
> )
> a1 += df.dot(v_mt, psi_H) * dx_heart
> 
> L1 = -I_ion * psi_H * dx_heart + v_mtm1 * psi_H * dx_heart
> 
> # eq2
> a2 = (
> M_i * df.inner(df.grad(v_mt), df.grad(psi_HaT)) * dx_heart
> + (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_heart
> + M_t * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_all(0)
> )
> 
> L2 = df.Constant(0) * psi_HaT * dx_all
> 
> a = a1 + a2
> L = L1 + L2
> 
> # some Boundary condition to test integration (phi otherwise wit non-zero nullspace)
> u_D = df.Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)
> 
> def boundary(x, on_boundary):
> return on_boundary
> 
> bc = df.DirichletBC(HaT, u_D, boundary)
> 
> sol = df.Function(W)
> df.solve(a == L, sol, bc)
> 
> ```

> **Working example**
>
> ```python
> import dolfin as df
> 
> mesh = df.UnitSquareMesh(20, 20)
> interior = df.CompiledSubDomain("std::max(fabs(x[0] - 0.5), fabs(x[1] - 0.5)) < 0.25")
> subdomains = df.MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
> # mark domain
> for cell in df.cells(mesh):
> x = cell.midpoint().array()
> subdomains[cell] = int(interior.inside(x, False))
> 
> # create interior mesh
> mesh_heart = df.MeshView.create(subdomains, 1)
> 
> # Function spaces
> H = df.FunctionSpace(mesh_heart, "Lagrange", 1)
> HaT = df.FunctionSpace(mesh, "Lagrange", 1)
> W = df.MixedFunctionSpace(H, HaT)
> 
> # Variational form
> dx_heart = df.Measure("dx", mesh_heart)
> dx_all = df.Measure("dx", mesh, subdomain_data=subdomains)
> 
> psi_H, psi_HaT = df.TestFunctions(W)
> v_mt, phi_et = df.TrialFunctions(W)
> 
> # last timestep v and ionic current
> # can this be done like this?
> v_mtm1 = df.Function(H)
> I_ion = df.Function(H)
> 
> # conductivities
> M_i = df.Constant(1)
> M_e = df.Constant(1)
> M_t = df.Constant(1)
> 
> # eq1
> a1 = (
> M_i * df.inner(df.grad(v_mt), df.grad(psi_H)) * dx_heart
> + M_i * df.inner(df.grad(phi_et), df.grad(psi_H)) * dx_heart
> )
> a1 += df.dot(v_mt, psi_H) * dx_heart
> 
> L1 = -I_ion * psi_H * dx_heart + v_mtm1 * psi_H * dx_heart
> 
> # eq2
> a2 = (
> M_i * df.inner(df.grad(v_mt), df.grad(psi_HaT)) * dx_heart
> + (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_all(1)
> + M_t * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_all(0)
> )
> 
> L2 = df.Constant(0) * psi_HaT * dx_all
> 
> a = a1 + a2
> L = L1 + L2
> 
> # some Boundary condition to test integration (phi otherwise wit non-zero nullspace)
> u_D = df.Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)
> 
> def boundary(x, on_boundary):
> return on_boundary
> 
> bc = df.DirichletBC(HaT, u_D, boundary)
> 
> sol = df.Function(W)
> df.solve(a == L, sol, bc)
> 
> ```

Again, Cécile, thank you for answering and for your efforts 🙂

Baltasar

---

_[View the full topic](https://fenicsproject.discourse.group/t/unable-to-successfully-call-petsc-function-matsetvalueslocal/8172)._