Anders Logg wrote:
> On Mon, Feb 18, 2008 at 02:54:38PM +0100, Kristen Kaasbjerg wrote:
>
>> Anders Logg wrote:
>>
>>> On Mon, Feb 18, 2008 at 02:29:23PM +0100, Kristen Kaasbjerg wrote:
>>>
>>>
>>>> Anders Logg wrote:
>>>>
>>>>
>>>>> On Mon, Feb 18, 2008 at 02:15:47PM +0100, Kristen Kaasbjerg wrote:
>>>>>
>>>>>
>>>>>
>>>>>> Anders Logg wrote:
>>>>>>
>>>>>>
>>>>>>
>>>>>>> Nice. The obvious thing would be to implement this in DiscreteFunction
>>>>>>> and map that to the function call
>>>>>>>
>>>>>>> virtual void Function::eval(real* values, const real* x) const;
>>>>>>>
>>>>>>> so that any Function (discrete, constant or user-defined) can be
>>>>>>> evaluated at an arbitrary point.
>>>>>>>
>>>>>>> It should be possible to implement this for any kind of element, and
>>>>>>> the code will look about the same as the code you have done for simple
>>>>>>> elements.
>>>>>>>
>>>>>>> We might add some kind of caching so that evaluation at multiple
>>>>>>> points that lie close to each other is efficient. (But maybe GTS is
>>>>>>> smart and handles this already.)
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>> ok guys, I have made a dirty hack in the C++ Function class
>>>>>> in order to get the desired functionality. Looks very much like
>>>>>> Dags code. Could I ask you to take a quick look at it (see below)
>>>>>> to see if I have done anything alarming. So now both the cell searching
>>>>>> and the function evaluation can be done from python (and perhaps be
>>>>>> condensed into one function call if desired) and it
>>>>>> seems to work.
>>>>>> Thanks for your help along the way.
>>>>>> Kristen
>>>>>>
>>>>>> -----------------------------------------------------------------------------------------------------------
>>>>>> void Function::my_eval(real* values, const real* x,
>>>>>> const ufc::cell& ufc_cell,
>>>>>> const ufc::finite_element& finite_element,
>>>>>> Cell& cell)
>>>>>> {
>>>>>> if (!f)
>>>>>> error("Function contains no data.");
>>>>>> //step #1: get expansion coefficient on the cell
>>>>>> uint n = finite_element.space_dimension();
>>>>>> real* coefficients = new real[n];
>>>>>> this->interpolate(coefficients,ufc_cell,finite_element,cell);
>>>>>>
>>>>>> //step #2: multiply with basis functions on the cell
>>>>>> real* basis_val = new real[finite_element.value_dimension(0)];
>>>>>> for(uint i=0; i<n; i++)
>>>>>> {
>>>>>> finite_element.evaluate_basis(i,basis_val,x,ufc_cell);
>>>>>> values[0] += basis_val[0]*coefficients[i];
>>>>>> }
>>>>>> }
>>>>>>
>>>>>>
>>>>>>
>>>>> Looks about right, but remember to delete the pointers coefficients
>>>>> and basis_val.
>>>>>
>>>>> Extending this to non-simple elements should be fairly simple. Add
>>>>> something like this:
>>>>>
>>>>> // Compute size of value (number of entries in tensor value)
>>>>> uint size = 1;
>>>>> for (uint i = 0; i < finite_element->value_rank(); i++)
>>>>> size *= finite_element->value_dimension(i);
>>>>>
>>>>> Then iterate over the number of values for each basis function (size),
>>>>> not only the first.
>>>>>
>>>>> Then we just need to include finding the element (using
>>>>> IntersectionDetector) in this function, remove the arguments ufc_cell,
>>>>> finite_element and cell, and then this can be added to DiscreteFunction.
>>>>>
>>>>>
>>>>>
>>>>>
>>>> Is the ufc::finite_element available in the Function class ?
>>>> Kristen
>>>>
>>>>
>>> No, but it's available in DiscreteFunction.
>>>
>>>
>>>
>> Ok, should I try to implement as much of this function as I can ?
>>
>
> Yes, that would be nice.
>
>
>> How are the Function and DiscreteFunction classes related and
>> what type is the FEM solution you get out from dolfin ?
>>
>
> It's a so-called envelope-letter design (with a twist).
>
> Basically, Function acts as the front-end for users, but does
> everything internally by calls to a pointer to a GenericFunction.
> This pointer is instantiated to either a DiscreteFunction,
> UserFunction or ConstantFunction depending on the arguments to the
> constructor of Function.
>
> So when you call u.eval() for a Function, then you call
> Function::eval(), which in turn calls GenericFunction::eval(), which
> is overloaded by for example DiscreteFunction::eval() depending on the
> representation of the function.
>
> In Function, you need to do something like
>
> void Function::eval(real* values, const real* x) const
> {
> if (!f)
> error("Function contains no data.");
> f->eval(values, x);
> }
>
> Then add eval() to the GenericFunction interface and implement eval()
> in DiscreteFunction, UserFunction (should return the same error as in
> Function now...) and ConstantFunction.
>
> See if you can find your way around...
>
>
Ok, have it implemented for simple elements now.
Seems to work for all subclasses of GenericFunction.
Should I bundle what I have and send it to you ?
I'm a little uncertain on the following things:
- point on boundary between two or more elements - is it then enough to
do the evaluation in one of the elements ?
- for UserFunction eval (which is only called when not overloaded) calls
the eval(x) of class Function (via f->eval(x)) to see if that then has
been overloaded, correct ?
- In DiscreteFunction.eval: a smart way to create the point (needed in
gts) depending of the spacial dimension of the finite_element.
Kristen
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev
