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;
}
------------------------------------------------------------------------------------------------------------
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev
