# Petsc error with calling dof\_to\_vertex

**URL:** https://fenicsproject.discourse.group/t/petsc-error-with-calling-dof-to-vertex/6237
**Category:** Uncategorized
**Created:** [July 15, 2021, 4:22pm UTC](https://fenicsproject.discourse.group/t/petsc-error-with-calling-dof-to-vertex/6237 "2021-07-15T16:22:52Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![liuyun](https://avatars.discourse-cdn.com/v4/letter/l/a8b319/32.png) [@liuyun](https://fenicsproject.discourse.group/u/liuyun)
#### Post date: [July 15, 2021, 4:22pm UTC](https://fenicsproject.discourse.group/t/petsc-error-with-calling-dof-to-vertex/6237/1 "2021-07-15T16:22:52Z")

</div>

I am trying to load a set of data using numpy and converting it to fenics mesh before solving the poisson equation. However, the following error is encountered

> Traceback (most recent call last):  
> File “sample\_poisson.py”, line 92, in   
> coeff0.vector().set\_local(coeff0\_val[d2v])  
> RuntimeError:
> 
> \*\*\* Error: Unable to successfully call PETSc function ‘VecSetValuesLocal’.  
> \*\*\* Reason: PETSc error code is: 63 (Argument out of range).  
> \*\*\* Where: This error was encountered inside /home/yl398/aur/fenics/dolfin/src/dolfin-2019.1.0.post0/dolfin/la/PETScVector.cpp.  
> \*\*\* Process: 0
> 
> * * *
> 
> \*\*\* DOLFIN version: 2019.1.0  
> \*\*\* Git changeset: acf365fd0a6e5ea80243a27621f7d9678f01d99f

I have attached a code here. The m\_h data contains 13^3=2197 entries, which matches the dimension of the mesh generated for the cube. There should be no problems with the boundary conditions. Any help is greatly appreciated.

```auto
from fenics import *
import numpy as np
from math import *

deg = 1
a = 12
## Create mesh and define function space
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(a, a, a), a, a, a)

m_h_input = np.genfromtxt('m_h.dat')

class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two slave edges
        return bool ((near(x[0], 0) or near(x[1], 0) or near(x[2], 0)) and 
            (not ((near(x[0], a) and near(x[2], a)) or 
                  (near(x[0], a) and near(x[1], a)) or
                  (near(x[1], a) and near(x[2], a)))) and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
    	#### define mapping for a single point in the box, such that 3 mappings are required
        if near(x[0], a) and near(x[1], a) and near(x[2],a):
            y[0] = x[0] - a
            y[1] = x[1] - a
            y[2] = x[2] - a
        ##### define mapping for edges in the box, such that mapping in 2 Cartesian coordinates are required
        if near(x[0], a) and near(x[2], a):
            y[0] = x[0] - a
            y[1] = x[1] 
            y[2] = x[2] - a      
        elif near(x[1], a) and near(x[2], a):
            y[0] = x[0] 
            y[1] = x[1] - a
            y[2] = x[2] - a
        elif near(x[0], a) and near(x[1], a):
            y[0] = x[0] - a
            y[1] = x[1] - a
            y[2] = x[2]         
        #### right maps to left: left/right is defined as the x-direction
        elif near(x[0], a):
            y[0] = x[0] - a
            y[1] = x[1]
            y[2] = x[2]
        ### back maps to front: front/back is defined as the y-direction    
        elif near(x[1], a):
            y[0] = x[0]
            y[1] = x[1] - a
            y[2] = x[2] 
        #### top maps to bottom: top/bottom is defined as the z-direction        
        elif near(x[2], a):
            y[0] = x[0]
            y[1] = x[1]
            y[2] = x[2] - a

## Define boundary condition
pbc = PeriodicBoundary()

V = FunctionSpace(mesh, 'Lagrange', deg, constrained_domain=pbc)
u = TrialFunction(V)
v = TestFunction(V)

d2v = dof_to_vertex_map(V)
coeff0_val = m_h_input
coeff0 = Function(V)
coeff0.vector().set_local(coeff0_val[d2v])

a = coeff0*inner(grad(u), grad(v))*dx
f = Constant(-6.0)
L = f*v*dx

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

```

---

<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: [July 17, 2021, 11:13am UTC](https://fenicsproject.discourse.group/t/petsc-error-with-calling-dof-to-vertex/6237/2 "2021-07-17T11:13:51Z")

</div>

As you are constraining your function space with a periodic boundary condition, you reduce the number of degrees of freedom in your space. Following your example, you can see this by printing the dim of a constrained and unconstrained function space:

```python
from fenics import *
import numpy as np
from math import *

deg = 1
a = 12
## Create mesh and define function space
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(a, a, a), a, a, a)

m_h_input = np.arange(2197, dtype=float)
#np.genfromtxt('m_h.dat')

class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two slave edges
        return bool ((near(x[0], 0) or near(x[1], 0) or near(x[2], 0)) and 
            (not ((near(x[0], a) and near(x[2], a)) or 
                  (near(x[0], a) and near(x[1], a)) or
                  (near(x[1], a) and near(x[2], a)))) and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
    	#### define mapping for a single point in the box, such that 3 mappings are required
        if near(x[0], a) and near(x[1], a) and near(x[2],a):
            y[0] = x[0] - a
            y[1] = x[1] - a
            y[2] = x[2] - a
        ##### define mapping for edges in the box, such that mapping in 2 Cartesian coordinates are required
        if near(x[0], a) and near(x[2], a):
            y[0] = x[0] - a
            y[1] = x[1] 
            y[2] = x[2] - a      
        elif near(x[1], a) and near(x[2], a):
            y[0] = x[0] 
            y[1] = x[1] - a
            y[2] = x[2] - a
        elif near(x[0], a) and near(x[1], a):
            y[0] = x[0] - a
            y[1] = x[1] - a
            y[2] = x[2]         
        #### right maps to left: left/right is defined as the x-direction
        elif near(x[0], a):
            y[0] = x[0] - a
            y[1] = x[1]
            y[2] = x[2]
        ### back maps to front: front/back is defined as the y-direction    
        elif near(x[1], a):
            y[0] = x[0]
            y[1] = x[1] - a
            y[2] = x[2] 
        #### top maps to bottom: top/bottom is defined as the z-direction        
        elif near(x[2], a):
            y[0] = x[0]
            y[1] = x[1]
            y[2] = x[2] - a

## Define boundary condition
pbc = PeriodicBoundary()

V_unc = FunctionSpace(mesh, 'Lagrange', deg)
V = FunctionSpace(mesh, 'Lagrange', deg, constrained_domain=pbc)
print(V_unc.dim(), V.dim())

```

which returns

```bash
2197 1797

```

---

<div class="post-metadata">

### Author: ![liuyun](https://avatars.discourse-cdn.com/v4/letter/l/a8b319/32.png) [@liuyun](https://fenicsproject.discourse.group/u/liuyun)
#### Post date: [July 19, 2021, 10:44am UTC](https://fenicsproject.discourse.group/t/petsc-error-with-calling-dof-to-vertex/6237/3 "2021-07-19T10:44:58Z")

</div>

Thank you very much for the explanation. That definitely helps a lot. I am having difficulties determining the exact degree of freedom using the periodic boundary conditions.

Initially there were a^3 degrees of freedom, but I’m mapping the 3 faces of the cube to another 3 faces, reducing the degree of freedom by 3_a^2. I’m mapping the edges separately which sort of avoids the double counting, so adding an addition 6_a degrees of freedom (since there are 12 edges).

So that should be a^3-3_a^2+6_a = 1768. This is different from the 1797 returned by the code. What am I doing wrong? Is the periodic boundary condition defined correctly?

---

<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: [July 19, 2021, 12:58pm UTC](https://fenicsproject.discourse.group/t/petsc-error-with-calling-dof-to-vertex/6237/4 "2021-07-19T12:58:53Z")

</div>

I would suggest you start with a slightly simpler mapping, (for instance map one facet to another facet), and try to understand what is going on from that mapping. I do not have time to double check that all your mappings are correct.

However, I think there might be a bug in DOLFIN when it comes to periodic mappings, as the number of dofs are not as expected in this MWE:

```python
from fenics import *
deg = 1
a = 2
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(a, a, a), a, a, a)

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return near(x[0], 0, 1e-6)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - 1.0
        y[1] = x[1]
        y[2] = x[2]

## Define boundary condition
pbc = PeriodicBoundary()

V_unc = FunctionSpace(mesh, 'Lagrange', deg)
V = FunctionSpace(mesh, 'Lagrange', deg, constrained_domain=pbc)

print(V_unc.dim(), V.dim())
print(V.tabulate_dof_coordinates())
print(V.dim(),(a+1)**3-(a+1)**2, (a+1)**3, V_unc.dim())
assert(V.dim() == (a+1)**3-(a+1)**2)

```

One thing you can always do to backward engineer the location of the input data is to tabulate dof coordinates (works on CG-1 spaces), and find the matching vertex.

You could also try out my implementation of periodic boundary conditions in dolfinx, namely:

> <https://github.com/jorgensd/dolfinx_mpc/blob/master/python/demos/demo_periodic_geometrical.py>

---

<div class="post-metadata">

### Author: ![liuyun](https://avatars.discourse-cdn.com/v4/letter/l/a8b319/32.png) [@liuyun](https://fenicsproject.discourse.group/u/liuyun)
#### Post date: [July 26, 2021, 12:58pm UTC](https://fenicsproject.discourse.group/t/petsc-error-with-calling-dof-to-vertex/6237/5 "2021-07-26T12:58:09Z")

</div>

Thank you very much for the suggestions. I will look at the example carefully.

It seems that I would need to install dolfinx for to use the code snippets? Is dolfinx stable for daily use?

---

<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: [July 26, 2021, 1:02pm UTC](https://fenicsproject.discourse.group/t/petsc-error-with-calling-dof-to-vertex/6237/6 "2021-07-26T13:02:27Z")

</div>

The first example is using dolfin, not dolfinx.

My implementation in dolfinx of periodic conditions are regulary updated to use the latest syntax of dolfinx. We have one release of dolfinx,

> **[Releases · FEniCS/dolfinx](https://github.com/FEniCS/dolfinx/releases)**
>
> Next generation FEniCS problem solving environment - FEniCS/dolfinx

Which is matched by:

> **[Releases · jorgensd/dolfinx\_mpc](https://github.com/jorgensd/dolfinx_mpc/releases)**
>
> Extension for dolfinx to handle multi-point constraints. - jorgensd/dolfinx\_mpc
