FEniCS Project

Using dolfin_function2BoxField to plot a solution on a structured grid from an solution on an imported mesh

Asked by caleb on 2010-04-21

Hello,

I have written a simple program to solve the Laplace equation over a triangular domain. The triangle has vertices at (1,1), (1,10), (10,1). I have created a mesh with Gmsh using a characteristic length of 1/5. My question is how to I use dolfin_function2BoxField to project my solution on a structured grid because I really want to use easyviz. The xml mesh is below:

<?xml version="1.0" encoding="UTF-8"?>

The script that I have written is as followed

import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import pylab as pl
import os
import scitools.easyviz as ev
from xml.etree import ElementTree as ET
from dolfin import *
from scitools.BoxField import *

"""
PHYS105: Homework 6
Problem 2

Solve Laplace's equation to calculate the steady-state
temperature distribution of a copper plate
"""

### Define problem parameters
lx = 2 # length of side x
ly = 1 # length of side y
d = 2 # degree of finite element, Lagrange polynominal

### Define boundary conditions
# Sub domaion for insulated boundary
# (mark whole boundary, UL and UB will overwrite)

class Noslip(SubDomain):
def inside(self, x, on_boundary):
return on_boundary

# Sub domain for upper left un-insolated section
class UL(SubDomain):
def inside(self, x, on_boundary):
return (abs(x[0] - 1.0) < DOLFIN_EPS and abs(x[1] - 5.0) > DOLFIN_EPS) and on_boundary

# Sub domain for bottom un-insolated section
class UB(SubDomain):
def inside(self, x, on_boundary):
return abs(x[1] - 1.0) < DOLFIN_EPS and on_boundary

### Read mesh and create function space
mesh = Mesh("triangle-test2.xml.gz")
V = FunctionSpace(mesh, 'CG', d)

### Mark subdomains
# Create mesh function over the cell facets
sub_domains = MeshFunction("uint", mesh, mesh.topology().dim() - 1)

# Mark all facets as sub domain 2
sub_domains.set_all(2)

# Mark upper-left Dirchelet boundary as sub domain 0
ul = UL()
ul.mark(sub_domains, 0)

# Mark bottom Dirchelet boundary as sub domain 1
ub = UB()
ub.mark(sub_domains, 1)

# Define boundary functions
left_DBC = Constant(200) # Upper left Dirchelet BC
bottom_DBC = Constant(100) # Bottom Dirchelet BC
g = Constant(0) # Everywhere else Neumann BC

bc = [DirichletBC(V, left_DBC, sub_domains, 0),
DirichletBC(V, bottom_DBC, sub_domains, 1)]

### Define variational problem
v = TestFunction(V)
u = TrialFunction(V)
f = Constant(0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx - g*v*ds(3)

### Compute solution
A = assemble(a, exterior_facet_domains=sub_domains)
b = assemble(L, exterior_facet_domains=sub_domains)
for condition in bc: condition.apply(A, b)

u = Function(V)
solve(A, u.vector(), b, 'lu')

### Project solution onto a mesh using one dimmensional elements
if u.ufl_element().degree() == 1:
u2 = u
else:
u2 = interpolate(u, FunctionSpace(mesh,'CG',1))

### Plot solution on structured grid
u_box = dolfin_function2BoxField(u2, mesh, (10,10), uniform_mesh=True)

ev.figure()
ev.mesh(u_box.grid.coorv[X], u_box.grid.coorv[Y], u_box.values,
title='mesh plot of U')

I can see that my problem comes from the divisions tuple because I have been essentially guessing how many divisions I would have using linspace(1,.2,10+.2) which might not have been the correct approach at all. How could the number of divisions be ascertained given just the mesh file as an input. The output of the script is

python problem2-test1.py
scitools.easyviz backend is gnuplot
Matrix of size 4610 x 4610 has 51872 nonzero entries.
Applying boundary conditions to linear system.
Applying boundary conditions to linear system.
Solving linear system of size 4610 x 4610 (PETSc LU solver, umfpack).
Traceback (most recent call last):
  File "problem2-test1.py", line 94, in <module>
    u_box = dolfin_function2BoxField(u2, mesh, (10,10), uniform_mesh=True)
  File "/usr/lib/pymodules/python2.6/scitools/BoxField.py", line 280, in dolfin_function2BoxField
    raise ValueError('Vector field (nodal_values) has length %d, there are %d grid points, and this does not match with %d components' % (nodal_values.size, grid.npoints, ncomponents))
ValueError: Vector field (nodal_values) has length 1191, there are 121 grid points, and this does not match with 9 components

Basically I just want to plot my solution using easyviz without having to export the solution and re-read the xml file like I had been doing before I started using dolfin_function2BoxField.

Thank you

Question information

Language:: English Edit question

Status:: Solved

For:: FEniCS Project Edit question

Assignee:: No assignee Edit question

Solved by:: Anders Logg

Solved:: 2010-04-29

Last query:: 2010-04-29

Last reply:: 2010-04-24

Link existing bug

Revision history for this message

Hans Petter Langtangen (hpl) said on 2010-04-24:

Wed, 21 Apr caleb wrote:
> New question #108000 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/108000
>
> Hello,
>
> I have written a simple program to solve the Laplace equation over a triangular domain. The triangle has vertices at (1,1), (1,10), (10,1). I have created a mesh with Gmsh using a characteristic length of 1/5. My question is how to I use dolfin_function2BoxField to project my solution on a structured grid because I really want to use easyviz. The xml mesh is below:

I would not recommend to use dolfin_function2BoxField in this case -
just visualize the data on the triangular domain directly using Viper
or VTK directly. If you insist on working with a structured,
rectangular mesh you may create such a mesh and set essential boundary
conditions at all nodes in the inactive triangle.

I'm not sure if it would work to compute u on your triangular domain
and then create a rectangular mesh m, an associated
V2=FunctionSpace(m, 'CG', 1) and then do a project(u, V2). You have
to test it.

Hans Petter

Revision history for this message

Anders Logg (logg) said on 2010-04-24:

On Sat, Apr 24, 2010 at 07:27:35PM -0000, Hans Petter Langtangen wrote:
> Question #108000 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/108000
>
> Status: Open => Answered
>
> Hans Petter Langtangen proposed the following answer:
> Wed, 21 Apr caleb wrote:
> > New question #108000 on FEniCS Project:
> > https://answers.launchpad.net/fenics/+question/108000
> >
> > Hello,
> >
> > I have written a simple program to solve the Laplace equation over a triangular domain. The triangle has vertices at (1,1), (1,10), (10,1). I have created a mesh with Gmsh using a characteristic length of 1/5. My question is how to I use dolfin_function2BoxField to project my solution on a structured grid because I really want to use easyviz. The xml mesh is below:
>
> I would not recommend to use dolfin_function2BoxField in this case -
> just visualize the data on the triangular domain directly using Viper
> or VTK directly. If you insist on working with a structured,
> rectangular mesh you may create such a mesh and set essential boundary
> conditions at all nodes in the inactive triangle.
>
> I'm not sure if it would work to compute u on your triangular domain
> and then create a rectangular mesh m, an associated
> V2=FunctionSpace(m, 'CG', 1) and then do a project(u, V2). You have
> to test it.
>
> Hans Petter

It is possible to just interpolate the function directly to a
rectangular mesh, if you set the option "allow_extrapolation" which
will extrapolate values in points outside domain instead of giving an
error message that the point was not found. Try this:

from dolfin import *

parameters["allow_extrapolation"] = True

circle = UnitCircle(16)
square = UnitSquare(16, 16)

Vc = FunctionSpace(circle, "CG", 1)
Vs = FunctionSpace(square, "CG", 1)

vc = project(Expression("sin(2*x[0])"), Vc)
vs = interpolate(vc, Vs)

plot(vc)
plot(vs)

interactive()

--
Anders

Revision history for this message

caleb (caselim) said on 2010-04-29:

Thanks Anders Logg, that solved my question.

To post a message you must log in.

Ask a question

Edit question

FEniCS Project

Using dolfin_function2BoxField to plot a solution on a structured grid from an solution on an imported mesh

Question information

Related bugs

Related FAQ:

Subscribers