# Assigning nodal values to a function only at the boundary nodes

I have been struggling about having a list of node indices:

nodes = [1, 10, 22, ..., 32]

of nodes at the boundary and an associated list of values

node_values = [10.2, 1.2,0.1, ..., -0.2]

and how to assign this to a function defined on the mesh, in order to use it as boundary conditions, instead of using an analytical function or a constant.

After searchin thoroughly I manage to find a solution that works for low order basis functions, for both scalars and vector. I put here my solution to see if this is the way to do it and how can I extend it to higher order basis functions. Also, if it helps someone:

import numpy

from time import time

from dolfin import *

def dofs_to_

"""

Compute the mapping from vertices to degrees of freedom only at the boundary.

This saves time, since we do not have to loop over all the cells of the domain.

If the for loop was replaced by a simple c++ function call, it would be much faster.

INPUTS ::

mesh :: The mesh where we want to compute the mapping between vertices

and degrees of freedom.

"""

# get the connectivity of the mesh

mesh.init(1) # initialize edges

mesh_topology = mesh.topology() # get the mesh topology

conn12 = mesh_topology(1,2) # get the cells to which an edge belongs

# get the number of edges

n_edges = boundaryFunctio

# get the indices of the boundary edges, as given by the boundaryFunction

boundaryEdges = numpy.arange(

boundaryEdges = boundaryEdges[

# Define the function spaces associated to the mesh

# for now only CG of order 1 works for sure

V = FunctionSpace(mesh, "CG", 1)

Vv = VectorFunctionS

# Allocate memory space of the array that stores the mapping from

# vertex number to degree of freedom number for both scalar and vector

# valued functions (note that for vector valued functions we need to indices

# since the x component and the y component have different indices)

vert_to_dofs = numpy.zeros(

vectordofs_

# Get the degree of freedom maps

dm = V.dofmap()

dms = [Vv.sub(i).dofmap() for i in range(2)] # the i=0 refers to x component

# get the cells of the mesh

mesh_cells = mesh.cells()

# Since we are only interested in the boundary vertices we just loop over

# the boundary edges as given by the MeshFunction and stored in boundaryEdges

# each time we loop we get the numbering of the degrees of freedom of the

# vertices of the cell that contains the boundary edge and we construct

# the mapping from vertex number to degree of freedom number

for edge in boundaryEdges:

# get the cell index to which the edge belongs

cell_ind = conn12(edge)[0]

# get the vertices of the cell

vert_inds = mesh_cells[

# add the degree of freedom numbering to the vertices

# now we do the same but for vector valued quantites

# since, in this case, our vector have two quantities (the x and y components)

# we have to iterate over the two sub dof maps

for i, (dms_i, dmcs_i) in enumerate(zip(dms, dms)):

# Return vertex to dof mapping

return vert_to_dofs, vectordofs_to_vert

# this is an example of how to assign values at the nodes for both a scalar

# and a vector valued function to a function

# create a unit square mesh (any other mesh works)

mesh = UnitSquareMesh(

# -------

# this part is optional, if you already know the indices of the vertices at

# the boundary just skip this part

# define the boundary mesh

boundaryMesh = BoundaryMesh(mesh)

# if you already don't have the boundary nodes indices

# get the boundary nodes indices

nodes = boundaryMesh.

# -------

# -------

# This part is just to generate a boundary function which assigns a value

# boundary_value to the boundary edges

boundary_value = 10

# get the global number of the boundary edges

edge_map = boundaryMesh.

# get the number of edges at the boundary

n_boundaryEdges = edge_map.size

# define the boundary function which defines the boundary

boundaryFunction = MeshFunction(

n_edges = boundaryFunctio

# allocate memory space to boundary values

boundary_values = numpy.zeros(

# assign the boundary_value to the array of boundary_values

boundary_

# assign the boundary_values to the boundary_function

boundaryFunctio

# -------

# generate some random values for the values to assign to the scalar function at the nodes

node_values_scalar = numpy.random.

# generate some random values for the values to assign to the vector function at the nodes

node_values_vector = numpy.random.

start = time()

# get the mapping of vertices at boundary to degrees of freedom

vert_to_dofs, vectordofs_to_vert = dofs_to_

print "It took :: %fs" %(time()-start)

# assign vertex values at the boundary for the scalar function

# generate the function space

V = FunctionSpace(

# generate the scalar function

myFunction = Function(V)

# get the vector containing the data values at the degrees of freedom

myFunction_vector = myFunction.vector()

# use the mapping from vertices (nodes) to degrees of freedom to assign the

# vertex values of the scalar function

myFunction_

plot(myFunction)

# assign vertex values at the boundary for the vector function

# generate the function space

Vv = VectorFunctionS

# generate the vector function

myVectorFunction = Function(Vv)

# get the vector containing the data values at the degrees of freedom

myVectorFunctio

# use the mapping from vertices (nodes) to degrees of freedom to assign the

# vertex values of the scalar function

myVectorFunctio

myVectorFunctio

plot(myVectorFu

interactive()

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Artur Palha

- Solved:
- 2013-03-21

- Last query:
- 2013-03-21

- Last reply:
- 2013-03-20

Artur Palha (artur-palha) said : | #1 |

One more question:

The for loop I have in the function dofs_to_verts is slow because the loop is made in python. This could be much faster if I could simlpy replace:

dm.cell_

by

list_of_cells = [2, 4, 5, ..., 22]

dm.cell_

Jan Blechta (blechta) said : | #2 |

If you think it can be helpful to someone else you can push it to demo dir after you debug it. Probably nobody would debug here so long code.

Jan

Artur Palha (artur-palha) said : | #3 |

Jan,

I have extended this code to higher order 'CG' elements. Do I have permissions to add things to the demo dir?

-artur palha

Jan Blechta (blechta) said : | #4 |

On Thu, 21 Mar 2013 13:51:10 -0000

Artur Palha <email address hidden> wrote:

> Question #224725 on DOLFIN changed:

> https:/

>

> Status: Answered => Solved

>

> Artur Palha confirmed that the question is solved:

> Jan,

>

> I have extended this code to higher order 'CG' elements. Do I have

> permissions to add things to the demo dir?

>

> -artur palha

>

Seems interesting. Actually you can't upload directly. You have to

create your branch and propose changes you make to be merged into

trunk. Follow instructions here:

http://

Jan

Jan Blechta (blechta) said : | #5 |

On Thu, 21 Mar 2013 15:52:27 +0100

Jan Blechta <email address hidden> wrote:

> On Thu, 21 Mar 2013 13:51:10 -0000

> Artur Palha <email address hidden> wrote:

> > Question #224725 on DOLFIN changed:

> > https:/

> >

> > Status: Answered => Solved

> >

> > Artur Palha confirmed that the question is solved:

> > Jan,

> >

> > I have extended this code to higher order 'CG' elements. Do I have

> > permissions to add things to the demo dir?

> >

> > -artur palha

> >

>

> Seems interesting. Actually you can't upload directly. You have to

> create your branch and propose changes you make to be merged into

> trunk. Follow instructions here:

> http://

>

> Jan

Don't do this now and save your time. FEniCS project is migrating from

launchapd to bitbucket at the moment.

Jan