ufl team mailing list archive

[Johan Hake <hake.dev@xxxxxxxxx>]

You can try generating a complex expression:

ccode = """
class Conditional: public Expression
{
public:
      double limit;
      Conditional():Expression()
      {
        limit = 5.0;
      }

      void eval(dolfin::Array<double>& values, const 
dolfin::Array<double>& x) const
      {
        values[0] = 10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 
0.02);
        values[0] = values[0]>limit ? limit : values[0];
      }
    };
"""
f = Expression(ccode)
f_projected = project(f,V)

print "the maximum value of f is: %g"% np.max(f_projected.vector().array())

Here you see that the projected maximal value is much closer to the 
limit. This might be a bug in UFL expressions?

Johan


On 07/12/2012 05:41 PM, Yaakoub El Khamra wrote:
> New question #202953 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/202953
>
> Greetings
> I have an inquiry regarding conditional expressions. I apologize in advance if this is more of a UFL issue than a Dolfin issue. I am trying to enforce limits on an expression (physically: species saturation). I am using something simple as follows:
>
> f_limited = conditional(ge(f, limit), limit, f)
>
> What I have noticed is that if I project f_limited for a quick plot, the conditional is not strictly honored. I am assuming this is due to the projection and nothing more, but I just want to check. The following is example code setup to replicate the issue, the limit is set to 5.0 but the maximum value of f_limited_projected turned out to be 5.03624. Is there a recommended method to inspect the limits (min/max) of an expression? Thank you very much in advance.
>
>
> from dolfin import *
> import numpy as np
>
> # Create mesh and define function space
> mesh = UnitSquare(32, 32)
> V = FunctionSpace(mesh, "CG", 3)
>
> f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
>
> f_projected = project(f,V)
>
> print "the maximum value of f is: %g"% np.max(f_projected.vector().array())
>
> print "establishing a limit of 5 on the value of f"
> limit = Constant(5.0)
> f_limited = conditional(ge(f, limit), limit, f)
>
> f_limited_projected = project(f_limited, V)
>
> print "the maximum value of f_limited is: %g"% np.max(f_limited_projected.vector().array())