I apologise for that, I missed a couple of details when I prepared a
minimum test case. The variable should have been grad_u not Efield.
Fixing that, the test case (in 2-D) does run successfully.
I've clarified that my problem comes from trying to solve in 1D rather
than 2D, starting from a UnitIntervalMesh instead of UnitSquareMesh.
I've double checked that the attached script gives the error
Vertex index (7) out of range [0, 7).
after identifying the refinement conditions
2 cells out of 4 marked for refinement (50.0%).
Drew
On Wed, 2015-07-08 at 11:42 +0200, Jan Blechta wrote:
> Attached script does not reproduce alleged problem. Instead it fails
> with message
>
> Traceback (most recent call last):
> File "demo_adaptive-poisson.py", line 74, in <module>
> gamma = abs(Efield.vector().array())
> NameError: name 'Efield' is not defined
>
> Could you fix it to reproduce the problem?
>
> Jan
>
>
> On Mon, 6 Jul 2015 11:41:20 +0000
> Drew Parsons <[email protected]> wrote:
>
> > Hi Fenics Folk,
> >
> > I'm not expert in Fenics/Dolfin, I'm evaluating to see if it can
> > perform better than a calculation I've currently got implemented in
> > octave (matlab).
> >
> > My problem in octave is that it fails to find the solution to
> > Poisson's equation when the electric field is extremely large. It
> > uses an adaptive spacing (1D mesh) but it complains that it can't
> > make
> > intervals any finer and reached the iteration limit. Perhaps
> > Fenics
> > can do a better job?
> >
> > So I'm testing your adaptive Poisson example,
> > /usr/share/dolfin/demo/undocumented/adaptive-poisson/
> > (maybe the documented auto-adaptive-poisson is better for me, but
> > anyway).
> >
> > The demo works fine. It uses "gamma" as the criterion for mesh
> > refinement, based on curvatures.
> >
> > I've replaced gamma in the example with the absolute value of the
> > gradient,
> > grad_u = project(grad(u), VectorFunctionSpace(mesh, "Lagrange",
> > 1)) gamma = abs(Efield.vector().array())
> >
> > The error estimate E (E=gamma*gamma) doesn't make sense then, but
> > mesh
> > refinement should still work, right?
> >
> > Mesh refinement does work twice, but in the third iteration it
> > fails
> > with the error:
> > *** Error: Unable to add cell using mesh editor.
> > *** Reason: Vertex index (13) out of range [0, 11).
> > *** Where: This error was encountered inside MeshEditor.cpp.
> >
> > The value of cell_markers going into refine() seems sensible, "3
> > cells
> > out of 7 marked for refinement (42.9%)." So I don't understand why
> > refine() would fail like this.
> >
> > So I wanted to check with you if this is a normal error or not. If
> > it
> > is normal than please tell me so and I'll keep studying the Fenics
> > system.
> >
> > I'm attaching my python file, which is identical to the one in the
> > adaptive-poisson example, apart from the change in definition of
> > gamma
> > (and a couple of print statements).
> >
> > I'm using DOLFIN version: 1.5.0 on Debian unstable.
> >
> > Thanks,
> >
> > Drew Parsons
> >
> """This demo program solves Poisson's equation
- div grad u(x, y) = f(x, y)
on the unit square with source f given by
f(x, y) = exp(-100(x^2 + y^2))
and homogeneous Dirichlet boundary conditions.
Note that we use a simplified error indicator, ignoring
edge (jump) terms and the size of the interpolation constant.
"""
# Copyright (C) 2008 Rolv Erlend Bredesen
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
#
# Modified by Anders Logg 2008-2011
from __future__ import print_function
from dolfin import *
from numpy import array, sqrt
from math import pow
from six.moves import xrange as range
TOL = 5e-4 # Error tolerance
REFINE_RATIO = 0.50 # Refine 50 % of the cells in each iteration
MAX_ITER = 20 # Maximal number of iterations
# Create initial mesh
mesh = UnitIntervalMesh(4)
source_str = "exp(-100.0*pow(x[0], 2))"
source = eval("lambda x: " + source_str)
# Adaptive algorithm
for level in range(MAX_ITER):
# Define variational problem
V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)
u = TrialFunction(V)
f = Expression(source_str)
a = dot(grad(v), grad(u))*dx
L = v*f*dx
# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, DomainBoundary())
# Compute solution
u = Function(V)
solve(a == L, u, bc)
# Compute error indicators
# h = array([c.diameter() for c in cells(mesh)])
# K = array([c.volume() for c in cells(mesh)])
# R = array([abs(source([c.midpoint().x(), c.midpoint().y()])) for c in cells(mesh)])
# gamma = h*R*sqrt(K)
grad_u = project(-grad(u), VectorFunctionSpace(mesh, "Lagrange", 1))
gamma = abs(grad_u.vector().array())
# Compute error estimate
E = sum([g*g for g in gamma])
E = sqrt(MPI.sum(mesh.mpi_comm(), E))
print("Level %d: E = %g (TOL = %g)" % (level, E, TOL))
# Check convergence
if E < TOL:
info("Success, solution converged after %d iterations" % level)
break
# Mark cells for refinement
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
gamma_0 = sorted(gamma, reverse=True)[int(len(gamma)*REFINE_RATIO)]
gamma_0 = MPI.max(mesh.mpi_comm(), gamma_0)
print("gamma_0=%g" % gamma_0)
print(gamma)
for c in cells(mesh):
cell_markers[c] = gamma[c.index()] > gamma_0
# Refine mesh
mesh = refine(mesh, cell_markers)
# Plot mesh
plot(mesh)
# Hold plot
interactive()
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
