Extrapolate to a higher order basis function (curve fit)

V1 = FunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 2)
v1 = Function(V1)
... populate v1.vector()
v2 = interpolate(v1, V2)

But since they are on the same mesh, the linear space is
a subspace of the quadratic space, and you will get the
same piecewise function, only using more dofs.

Martin

On 23 November 2011 20:20, Charles Cook
<email address hidden> wrote:
> New question #179778 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/179778
>
>
> I'm trying to extrapolate a linear function space onto a quadratic function space, where the linear function space is populated from data.  Each DOF in the linear function space is set along with the mesh coordinate locations from data.  Then a quadratic function space is created on the same mesh, with of course more DOFs.  The quadratic function space is then set by extrapolating the linear function space.  This is done by interpolating on the linear function space to find the necessary DOFs for the quadratic function space? (I believe this to be the case).
>
> Is there a way I can 'fit' the second order function space to my data?  IE, curve fit?
>
> My motivation is to use the FEM basis functions to interpolate the data points.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Thank you Martin,

This is exactly what I am looking for, however, where may I find interpolate in C++?

This appears to be what I need?

http://fenicsproject.org/documentation/dolfin/dev/cpp/programmers-reference/function/FunctionSpace.html#FunctionSpace::interpolate__GenericVectorR.GenericFunctionCRC

Also, I just caught your comment... So the interpolation level would still be first order, just using more DOFs?

I don't suppose there is a way to get a higher level of interpolation? This would be key for a data driven coefficient with higher than first order basis functions... Otherwise the solution would be contaminated from the data?

Thank you Anders. This would be great! I tried implementing this, however, I get an erroneous extrapolation.

My test code:

// Setup data vectors
Vector xs(n);
Vector ys(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]);
}

// Setup the mesh to interpolate on
// Note the count is the number of elements, not nodes
Interval mesh (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.
Interpolate1::FunctionSpace V (mesh);

// copy the data into the function space
Function fy (V, ys);

Interpolate2::FunctionSpace V2 (mesh);
Function fy2 (V);

fy2.extrapolate(fy);

std::cout << "O(1): " << fy(2.2)
  << "\tO(2): " << fy2(2.2)
  << "\n";

where

Interpolate1:
# First order (linear) lagrange elements (polynomials)
element = FiniteElement("Lagrange", interval, 1)

Interpolate2:
# Second order lagrange elements (polynomials)
element = FiniteElement("Lagrange", interval, 2)

The output is:

O(1): 4.84 O(2): 6.05

where I would expect the result to be 4.84 for both cases.

Why is it *always* once you post an example you see the problem? I put the wrong function space into fy2. The corrected code which works wonderfully is included bellow. Note that I test at x=2.0 which is off a node and converges to 2^2 correctly.

Thank you!!!

// Setup data vectors
Vector xs(n);
Vector ys(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]);
}

// Setup the mesh to interpolate on
// Note the count is the number of elements, not nodes
Interval mesh (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.
Interpolate1::FunctionSpace V (mesh);

// copy the data into the function space
Function fy (V, ys);

Interpolate2::FunctionSpace V2 (mesh);
Function fy2 (V2);

fy2.extrapolate(fy);

std::cout << "O(1): " << fy(2)
  << "\tO(2): " << fy2(2)
  << "\n";

Thanks Anders Logg, that solved my question.

On Wed, Nov 23, 2011 at 08:45:46PM -0000, Charles Cook wrote:
> Question #179778 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/179778
>
> Status: Open => Solved
>
> Charles Cook confirmed that the question is solved:
> Why is it *always* once you post an example you see the problem? I put
> the wrong function space into fy2. The corrected code which works
> wonderfully is included bellow. Note that I test at x=2.0 which is off
> a node and converges to 2^2 correctly.
>
> Thank you!!!
>
> // Setup data vectors
> Vector xs(n);
> Vector ys(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]);
> }
>
>
> // Setup the mesh to interpolate on
> // Note the count is the number of elements, not nodes
> Interval mesh (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.
> Interpolate1::FunctionSpace V (mesh);
>
> // copy the data into the function space
> Function fy (V, ys);
>
>
> Interpolate2::FunctionSpace V2 (mesh);
> Function fy2 (V2);
>
> fy2.extrapolate(fy);
>
> std::cout << "O(1): " << fy(2)
> << "\tO(2): " << fy2(2)
> << "\n";

Good to hear that it works!

--
Anders

Just an updated code listing for the new shared pointer usage.

 // 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]);
 }

 // 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));;

 Function fy2 (V2);

 fy2.extrapolate(fy);

 std::cout << "O(1): " << fy(2.12)
   << "\tO(2): " << fy2(2.12)
   << "\n";

Hello,

I've started working against the development repositories and am running into an odd problem now. The same methods used above now throw a dolfin_error.

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** https://answers.launchpad.net/dolfin
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to compute extrapolation.
*** Reason: Not enough degrees of freedom on local patch to build extrapolation.
*** Where: This error was encountered inside Extrapolation.cpp.
*** Process: 0
*** -------------------------------------------------------------------------

I put a break right before the error was thrown in Extrapolation.cpp,

  // Compute size of linear system
  dolfin_assert(W.element());
  const uint N = W.element()->space_dimension();
  const uint M = unique_dofs.size();

  // Check size of system
  if (M < N)
  {
    dolfin_error("Extrapolation.cpp",
                 "compute extrapolation",
                 "Not enough degrees of freedom on local patch to build extrapolation");
  }

M was valued at 3, using FiniteElement("Lagrange", interval, 4), no problem there

N was oddly optimized out -- but if my understanding is correct, since both are intervals, N should be 1 for 1-D, making M>N?

Thank you!

Charles

The space_dimension() of an element is the dimension of the local finite element space, in other words, the number of local basis functions or degrees of freedom. For 4th order Lagrange on an interval, that dimension would be 4+1 = 5 > 3.

Thank you for the response Marie. The intention of the comparison then makes sense.

Extrapolating a 1st order Lagrange on an interval to a 4th order Lagrange on an interval should be

M = 3 (DOF on source as pulled from memory)
N = 5 (DOF on destination)

if (M < N) //false

But this condition always evaluates true for me.

Could

const uint N = W.element()->space_dimension();

be an error? (checking space_dimension not DOF?)

You might want to check that logic again: with M = 3 and N = 5: M < N, so it is not that suprising that this condition indeed evaluates to true.

Wow - you're absolutely right. I was thinking DOF for 4th order was greater than 1st and got it crossed by expectation.

Taking a step back... perhaps I am using the method wrong, how would I extrapolate from a 1st order langrange interval on a 4th order lagrange interval?

I have been doing something like

Function fv1(V1 ys);

Function fv4(V4);

fv4.extrapolate(fv1);

Where V1 and V4 are the respective forms.