Kristen Kaasbjerg wrote:
> Kristen Kaasbjerg wrote:
>> 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
>>
>>
>
> Hhhhmmm, something seems to have changed between version 0.7.1 and
> 0.7.2. When calling the eval function of class DiscreteFunction,
> the call to finite_element->evaluate_basis returns strange small
> numbers like 5.81e-268. Has anything changed between these two versions ?
> My implementation:
> --------------------------------------------------------------------------------------------------------
>
>
> void DiscreteFunction::eval(real* values, const real* x)
> {
> uint n = finite_element->space_dimension();
> //step #1: locate the cell that contains x
> // check also if x has the correct dimension
> Point p(x[0],x[1]); // generalize to arbitrary dimension
> IntersectionDetector id(mesh);
> Array<uint> cells;
> id.overlap(p,cells);
> uint cell_index = cells[0];
> Cell cell(mesh,cell_index);
> UFCCell ufc_cell(cell);
>
> //step #2: get expansion coefficient on the cell
> real* coefficients = new real[n];
> this->interpolate(coefficients,ufc_cell,*finite_element);
>
> //step #3: multiply with basis functions on the cell
> real* basis_val = new real[finite_element->value_dimension(0)];
> real value(0.);
> for(uint i=0; i<n; i++)
> {
> finite_element->evaluate_basis(i,basis_val,x,ufc_cell);
> value += basis_val[0]*coefficients[i];
> cout << basis_val[0] << endl;
> }
> values[0] = value;
> delete [] coefficients;
> delete [] basis_val;
> }
> ------------------------------------------------------------------------------------------------------------
>
>
I think I'll need some expertise to sort out this strange behavior.
I think it might be related to the python wrappings. At least, when
I make a C++ version of my demo, the correct values are returned
from the evaluate_basis function call. Whereas in python undefined
behavior is observed.
Kristen
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev
