extracting solution values at nodes for arbitrary degree finite element

On Sun, Apr 18, 2010 at 07:57:51AM -0000, caleb wrote:
> New question #107724 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/107724
>
> Hello,
>
> After using Dolfin to solve a variational problem, there are many times where I want the solution at the mesh nodes for plotting on a structured grid for example. When I use a finite element of degree 1, I have no problem using the dolfin_function2BoxField(), but when I use a degree of higher order I find that if I assign u to contain my solution I obtain a u.vector().array() with a size of d times the number of grid points so dolfin_function2BoxField() throws an error message. Is there some option that I should pass to dolfin_function2BoxField() to inform it of the degree of my finite element?

What is dolfin_function2BoxField?

It's not part of DOLFIN.

--
Anders

Sun, 18 Apr Anders Logg wrote:
> Question #107724 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/107724
>
> Status: Open => Answered
>
> Anders Logg proposed the following answer:
> On Sun, Apr 18, 2010 at 07:57:51AM -0000, caleb wrote:
> > New question #107724 on FEniCS Project:
> > https://answers.launchpad.net/fenics/+question/107724
> >
> > Hello,
> >
> > After using Dolfin to solve a variational problem, there are many times where I want the solution at the mesh nodes for plotting on a structured grid for example. When I use a finite element of degree 1, I have no problem using the dolfin_function2BoxField(), but when I use a degree of higher order I find that if I assign u to contain my solution I obtain a u.vector().array() with a size of d times the number of grid points so dolfin_function2BoxField() throws an error message. Is there some option that I should pass to dolfin_function2BoxField() to inform it of the degree of my finite element?

Can you specify the error message and provide a piece of code that
triggers the error? Without it, it's difficult to solve the problem.

> What is dolfin_function2BoxField?
>
> It's not part of DOLFIN.

dolfin_function2BoxField is a function introduced in the FEniCS tutorial
to map DOLFIN finite element fields onto "finite difference" grids
for visualization and analysis, by a number of tools.

Hans Petter

On Sun, Apr 18, 2010 at 03:48:21PM -0000, Hans Petter Langtangen wrote:
> Question #107724 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/107724
>
> Hans Petter Langtangen proposed the following answer:
> Sun, 18 Apr Anders Logg wrote:
> > Question #107724 on FEniCS Project changed:
> > https://answers.launchpad.net/fenics/+question/107724
> >
> > Status: Open => Answered
> >
> > Anders Logg proposed the following answer:
> > On Sun, Apr 18, 2010 at 07:57:51AM -0000, caleb wrote:
> > > New question #107724 on FEniCS Project:
> > > https://answers.launchpad.net/fenics/+question/107724
> > >
> > > Hello,
> > >
> > > After using Dolfin to solve a variational problem, there are many times where I want the solution at the mesh nodes for plotting on a structured grid for example. When I use a finite element of degree 1, I have no problem using the dolfin_function2BoxField(), but when I use a degree of higher order I find that if I assign u to contain my solution I obtain a u.vector().array() with a size of d times the number of grid points so dolfin_function2BoxField() throws an error message. Is there some option that I should pass to dolfin_function2BoxField() to inform it of the degree of my finite element?
>
> Can you specify the error message and provide a piece of code that
> triggers the error? Without it, it's difficult to solve the problem.
>
> > What is dolfin_function2BoxField?
> >
> > It's not part of DOLFIN.
>
> dolfin_function2BoxField is a function introduced in the FEniCS tutorial
> to map DOLFIN finite element fields onto "finite difference" grids
> for visualization and analysis, by a number of tools.

Aha, it shows I haven't edited that part of the tutorial yet. ;-)

--
Anders

Hello,

The error output is as followed

Traceback (most recent call last):
  File "test1.py", line 98, in <module>
    u_box = dolfin_function2BoxField(u, mesh, (nx,ny), 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 2501, there are 651 grid points, and this does not match with 3 components

The Python script is below

import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import pylab as pl
import os
import mpl_toolkits.mplot3d.axes3d as axes3d
from xml.etree import ElementTree as ET
from dolfin import *
from scitools.BoxField import *
from scitools.easyviz import plot as curve_plot, \
     figure as new_figure, \
     contour, \
     surf as surf_plot, mesh as mesh_plot

# PHYS105: Test 1, question 1 (2D)
# Solve the diffusion equation on a hollow cylinder where there is no angular dependence on the heat flow

# Create mesh and define function space
nx = 20
ny = 30
d = 2
r1 = 0.5
r2 = 1.0
T = 10 # total simulation time
dt = 0.2 # time step

#mesh = UnitSquare(nx, ny)
mesh = Rectangle(r1,0,r2,1, nx, ny)
V = FunctionSpace(mesh, 'CG', d)

# Extract mesh coordinates
x = mesh.coordinates()[:,0]
y = mesh.coordinates()[:,1]

# Generate hollow circle, Theta radians,
# with inner radius r1 and outer radius r2
theta = pi/2

def cylinder(r, s):
    return [r*np.cos(theta*s), r*np.sin(theta*s)]

x_hat, y_hat = cylinder(x, y)
xy_hat_coor = np.array([x_hat, y_hat]).transpose()

mesh.coordinates()[:] = xy_hat_coor
X = mesh.coordinates()[:,0]
Y = mesh.coordinates()[:,1]

# Define boundary conditions
u_a = Expression('t', element = V.ufl_element(), degree =d)
u_b = Expression('100+40*t', element = V.ufl_element(), degree =d)

class DirichletBoundary0(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14 # tolerance for coordinate comparisons
        return on_boundary and ( (np.sqrt(x[0]**2 +x[1]**2) - r1) < tol)

class DirichletBoundary1(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14 # tolerance for coordinate comparisons
        return on_boundary and ( (np.sqrt(x[0]**2 +x[1]**2) - r2) < tol)

bc = [DirichletBC(V, u_b, DirichletBoundary1()), DirichletBC(V, u_a, DirichletBoundary0())]

# Define the initial condtion
u_a.t = 0
u_b.t = 0
u0 = Expression("200*(sqrt((x[0]*x[0] + x[1]*x[1]))-0.5)", element = V.ufl_element(), degree =d)
# Initial condition
u_prev = interpolate(u0, V)

# Define variational problem
v = TestFunction(V)
u = TrialFunction(V)
p = Expression('4.0', element = V.ufl_element(), degree =d)
f = Expression('0', element = V.ufl_element(), degree =d)
a = u*v*dx + p*dt*dot(grad(u), grad(v))*dx
L = (u_prev + dt*f)*v*dx

# Compute solution

t = dt
while t <= T:
    print 'time =', t
    u_a.t = t
    u_b.t = t
    problem = VariationalProblem(a, L, bc)
    u = problem.solve()
    t += dt
    u_prev.assign(u)

u_box = dolfin_function2BoxField(u, mesh, (nx,ny), uniform_mesh=True)
new_figure()
mesh_plot(u_box.grid.coorv[X], u_box.grid.coorv[Y], u_box.values,
          title='mesh plot of U')

plt.show()

Hello,

The problem seems to come from assigning V to have a degree > 1 if that helps at all

Mon, 19 Apr caleb wrote:
> Question #107724 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/107724
>
> caleb gave more information on the question:
> Hello,
>
> The problem seems to come from assigning V to have a degree > 1 if that
> helps at all

Can you provide a test program? It's very difficult for me to track
down any bug without a specific program example that fails and that I
can experiment with.

Hans Petter

Wonderful, I will try to incorporate the change this afternoon and hopefully mark this question as solved.

Thanks Hans Petter Langtangen, that solved my question.

hello,
 any one who can help me ,how to create mesh in irregular domain using Fenics.

Please start a new question.