Errors in Function.extrapolate with Interval Lagrange elements

Also, I am using the latest Ubuntu PPA binaries.

Note that if I do the following instead the extrapolations is stable, but only maintains the second order accuracy:

 fy2.extrapolate(fy);
 fy3.extrapolate(fy2);
 fy4.extrapolate(fy3);

It appears increasing the order more than one degree causes a problem?

The problem is that there aren't enough degrees of freedom on the
local patch (in 1D consisting of 3 intervals) to extrapolate from
degree 1 to degree 4. There are just enough degrees of freedom
to extrapolate from degree 1 to degree 3.

I have added a test for this now that will result in an error message
when there are too few degrees of freedom. This will be pushed to the
main branch as soon as the tests pass.

--
Anders

On Wed, Nov 30, 2011 at 12:45:49AM -0000, Charles Cook wrote:
> New question #180447 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/180447
>
> I've been testing extrapolate for interpolation purposes and noticed what appears to be a bug (code below). When extrapolating more than one order higher in interpolation levels (say Lagrange order 1 to Lagrange order 3) oscillations appear in the extrapolated higher order function.
>
> The following is a test code which generates data for y=x^5 at non-regular intervals, x_i = i*(1+*0.1), and tests extrapolating the data to higher basis functions at x=3. The following is the output:
>
>
> 3^5 : 243 <--- expected
> O(1): 298.676 O(2): 228 O(3): 234.062 O(4): 88.9269
>
> Note that the values start to diverge. The errors behave differently with different data and testing points, but this demonstrates the problem.
>
> The UFL files look like:
>
> element = FiniteElement("Lagrange", interval, 1) # Interpolate1.ufl
> element = FiniteElement("Lagrange", interval, 2) # Interpolate2.ufl
> element = FiniteElement("Lagrange", interval, 3) # Interpolate3.ufl
> element = FiniteElement("Lagrange", interval, 4) # Interpolate4.ufl
>
>
>
> Here's the code:
>
>
> // Setup data vectors
> int n = 9;
> // Setup data vectors
> boost::shared_ptr<Vector> xs(new Vector(n));
> boost::shared_ptr<Vector> ys(new Vector(n));
>
> // Populate test data
> for(int i=0; i<n; i++)
> {
> xs->setitem(i,i+i*0.1);
> ys->setitem(i,(*xs)[i]*(*xs)[i]*(*xs)[i]*(*xs)[i]*(*xs)[i]);
> }
>
> // Setup the mesh to interpolate on
> // Note the count is the number of elements, not nodes
> boost::shared_ptr<Mesh> mesh (new Interval(n-1,0,1));
> for(int i=0; i<n; i++)
> mesh->coordinates()[i] = (*xs)[i];
>
> // create the function space of order one.
> // this is done so that dof maps directly to
> // data points.
> boost::shared_ptr<Interpolate1::FunctionSpace> V
> (new Interpolate1::FunctionSpace(mesh));
>
> // copy the data into the function space
> Function fy (V, ys);
>
> boost::shared_ptr<Interpolate2::FunctionSpace> V2
> (new Interpolate2::FunctionSpace(mesh));;
> boost::shared_ptr<Interpolate3::FunctionSpace> V3
> (new Interpolate3::FunctionSpace(mesh));;
> boost::shared_ptr<Interpolate4::FunctionSpace> V4
> (new Interpolate4::FunctionSpace(mesh));;
>
> Function fy2 (V2);
> Function fy3 (V3);
> Function fy4 (V4);
>
> fy2.extrapolate(fy);
> fy3.extrapolate(fy);
> fy4.extrapolate(fy);
>
> std::cout << "3^5 : " << 3*3*3*3*3 << "\n"
> << "O(1): " << fy(3)
> << "\tO(2): " << fy2(3)
> << "\tO(3): " << fy3(3)
> << "\tO(4): " << fy4(3)
> << "\n";
>
>
> Thoughts?
>

Thank you!

Though, I note when testing with a degree 3 polynomial I also see oscillations. With the same method above using the following data:

data = [373.15 101300
    400 247000
    430 571000
    460 1172000];

I obtain the following values using a degree 3 polynomial:

intData3 = [374 76274.7
379 -18033.6
384 -33991
389 11707.5
394 102367
399 221292
404 278177
409 321704
414 370772
419 425895
424 487584
429 556352
434 304527
439 33056.5
444 -104500
449 -38021.3
454 302616
459 987532];

There are some dramatic oscilations in this?

Degree 2 looks better:

intData2 = [374 104055
379 122390
384 144372
389 169999
394 199271
399 232190
404 271205
409 316920
414 368994
419 427426
424 492216
429 563365
434 639805
439 724275
444 817257
449 918751
454 1.02876e+06
459 1.14727e+06];

Is there a way to directly set the DOFs for a higher order polynomial as is done for the order 1 polynomial? I'd like to assign them based on say a spline fit of the data.

This question was expired because it remained in the 'Open' state without activity for the last 15 days.