dolfin team mailing list archive

[Kristian Oelgaard <k.b.oelgaard@xxxxxxxxxx>]

Quoting Anders Logg <logg@xxxxxxxxx>:

> On Mon, Oct 05, 2009 at 11:18:19AM +0200, Kristian Oelgaard wrote:
> >
> > Hi,
> >
> > I'm trying to update my plasticity solver to use the latest version of
> DOLFIN
> > and I ran into a problem with the new design of Function/Expression etc.
> > I get the following error:
> >
> > terminate called after throwing an instance of 'std::runtime_error'
> >   what():  evaluate_basis() is not supported for QuadratureElement
> >
> > Program received signal SIGABRT, Aborted.
> > [Switching to Thread 0xb6a436d0 (LWP 21130)]
> > 0xb7f15430 in __kernel_vsyscall ()
> > (gdb) where
> > #0  0xb7f15430 in __kernel_vsyscall ()
> > #1  0xb6b726d0 in raise () from /lib/tls/i686/cmov/libc.so.6
> > #2  0xb6b74098 in abort () from /lib/tls/i686/cmov/libc.so.6
> > #3  0xb6d8f8f8 in __gnu_cxx::__verbose_terminate_handler () from
> > /usr/lib/libstdc++.so.6
> > #4  0xb6d8d7d5 in ?? () from /usr/lib/libstdc++.so.6
> > #5  0xb6d8d812 in std::terminate () from /usr/lib/libstdc++.so.6
> > #6  0xb6d8d94a in __cxa_throw () from /usr/lib/libstdc++.so.6
> > #7  0x080d6b39 in strain2d_0_finite_element_0::evaluate_basis
> (this=0x9ba9038,
> > i=0, values=0x9bad220, coordinates=0xbfb2f3c8, c=@0xbfb2f5d8)
> >     at
> >
>
/home/oelgaard/fenics_projects/fenics_apps/fenics-plasticity/demo/forms/2Dlinear_QE_forms/Strain2D.cpp:495
> > #8  0xb7cf2d55 in dolfin::Function::eval (this=0xbfb306b4,
> values=0xbfb2f3b0,
> > x=0xbfb2f3c8, ufc_cell=@0xbfb2f5d8, cell_index=0) at
> > ./dolfin/fem/FiniteElement.h:53
> > #9  0xb7cf9e77 in dolfin::Function::eval (this=0xbfb306b4,
> values=0xbfb2f3b0,
> > data=@0xbfb306fc) at dolfin/function/Function.cpp:275
> > #10 0xb7d00c77 in dolfin::GenericFunction::evaluate (this=0xbfb306f8,
> > values=0xbfb2f3b0, coordinates=0xbfb2f3c8, cell=@0xbfb2f5d8)
> >     at dolfin/function/GenericFunction.cpp:35
> > #11 0x080bc905 in plas2d_1_finite_element_2::evaluate_dof (this=0x9bb7558,
> i=0,
> > f=@0xbfb306f8, c=@0xbfb2f5d8)
> >     at
> >
>
/home/oelgaard/fenics_projects/fenics_apps/fenics-plasticity/demo/forms/2Dlinear_QE_forms/Plas2D.cpp:10433
> > #12 0xb7d00f4c in dolfin::GenericFunction::restrict_as_ufc_function
> > (this=0xbfb306f8, w=0x9bb7c38, element=@0x9bb7580, dolfin_cell=@0xbfb2f520,
> >     ufc_cell=@0xbfb2f5d8, local_facet=-1) at
> ./dolfin/fem/FiniteElement.h:61
> > #13 0xb7cf2f8e in dolfin::Function::restrict (this=0xbfb306b4, w=0x9bb7c38,
> > element=@0x9bb7580, dolfin_cell=@0xbfb2f520, ufc_cell=@0xbfb2f5d8,
> > local_facet=-1)
> >     at dolfin/function/Function.cpp:490
> > #14 0xb7ceb0ce in dolfin::UFC::update (this=0xbfb2f5b0, cell=@0xbfb2f520)
> at
> > dolfin/fem/UFC.cpp:212
> > #15 0xb7cb2e85 in dolfin::Assembler::assemble_cells (A=@0x9bbb838,
> > a=@0xbfb30b18, ufc=@0xbfb2f5b0, domains=0x0, values=0x0) at
> > dolfin/fem/Assembler.cpp:159
> > #16 0xb7cb3475 in dolfin::Assembler::assemble (A=@0x9bbb838, a=@0xbfb30b18,
> > cell_domains=0x0, exterior_facet_domains=0x0, interior_facet_domains=0x0,
> >     reset_sparsity=<value optimized out>, add_values=<value optimized out>)
> at
> > dolfin/fem/Assembler.cpp:115
> > #17 0xb7cb35e2 in dolfin::Assembler::assemble (A=@0x9bbb838, a=@0xbfb30b18,
> > reset_sparsity=<value optimized out>, add_values=<value optimized out>)
> >     at dolfin/fem/Assembler.cpp:50
> > #18 0xb7eef1f8 in dolfin::PlasticityProblem::F (this=0xbfb2f814,
> b=@0x9bbb838,
> > x=@0x9bb0e48) at src/PlasticityProblem.cpp:92
> > #19 0xb7df9baf in dolfin::NewtonSolver::solve (this=0xbfb308d4,
> > nonlinear_problem=@0xbfb2f814, x=@0x9bb0e48) at
> dolfin/nls/NewtonSolver.cpp:76
> > #20 0x08096370 in main () at main.cpp:184
> >
> > It's pretty obvious what the error is, but before you say that we should
> just
> > implement evaluate_basis() for QuadratureElement (which is alway 1, if it's
> > evaluate exactly at the quadrature point, otherwise it's undefined) I would
> > like to know why it is necessary to call evaluate on a function that
> contains a
> > computed result. It should be possible to just pass this function to the
> Form
> > which will handle the interpolation when integrating the form.
> >
> > Kristian
>
> I don't understand why evaluate_basis is called. The assembler
> restricts all coefficients to the local element which results in a
> call to evaluate_dof which results in a call to evaluate. Why is
> Strain2D.cpp calling evaluate_basis on line 495? It should just return
> the value.

evaluate() ends up with a call to Function::eval() which calls evaluate_basis()

The attached forms will reproduce the error.

Kristian




> --
> Anders
>


// Copyright (C) 2006-2009 Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
// First added:  2006-02-07
// Last changed: 2009-10-05
//
// 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) = 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)
//
// and boundary conditions given by
//
//     u(x, y) = 0        for x = 0 or x = 1
// du/dn(x, y) = sin(5*x) for y = 0 or y = 1

#include <dolfin.h>
#include "Quad.h"
#include "Project.h"

using namespace dolfin;

// Source term (right-hand side)
class Source : public Expression
{
public:

  Source() : Expression(2) {}

  void eval(double* values, const double* x) const
  {
    double dx = x[0] - 0.5;
    double dy = x[1] - 0.5;
    values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
  }

};

// Boundary flux (Neumann boundary condition)
class Flux : public Expression
{
public:

  Flux() : Expression(2) {}

  void eval(double* values, const double* x) const
  {
    values[0] = sin(5*x[0]);
  }

};

// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
  bool inside(const double* x, bool on_boundary) const
  {
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS;
  }
};

int main()
{
  // Create mesh and function space
  UnitSquare mesh(32, 32);
  Quad::FunctionSpace V(mesh);

  // Define boundary condition
  Constant u0(mesh, 0.0);
  DirichletBoundary boundary;
  DirichletBC bc(V, u0, boundary);

  // Define variational problem
  Quad::BilinearForm a(V, V);
  Quad::LinearForm L(V);
  Source f;
  Flux g;
  L.f = f;
//  L.g = g;

  // Compute solution
  VariationalProblem problem(a, L, bc);
  Function u(V);
  problem.solve(u);

  Project::FunctionSpace Vp(mesh);
  Project::BilinearForm ap(Vp, Vp);
  Project::LinearForm Lp(Vp);
  Lp.u0 = u;

  VariationalProblem problemp(ap, Lp);
  Function U(Vp);
  problemp.solve(U);

  // Save solution in VTK format
  File file("poisson.pvd");
  file << U;

  // Plot solution
//  plot(u);

  return 0;
}

Attachment: Project.ufl
Description: Binary data

Attachment: Quad.ufl
Description: Binary data

Attachment: SConstruct
Description: Binary data