dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Hi everyone!

As you may have noticed, I've been working on some new demos in
src/demo/pde. As a consequence, I have fixed a couple of issues with
the Function class that we've had previously. In particular, you no
longer have to initalize a Vector for a Function used as initial value
in time-steppig.

Here's an illustration from the new convection-diffusion demo:

// Linear system and solver
Matrix A;
Vector x, b;
GMRES solver;
                                                                                                    
// Create functions
ConvectionDiffusion::BilinearForm::TrialElement element;
Function u0 = 0.0;
Function u1(x, mesh, element);
                                                                                                    
// Create forms
ConvectionDiffusion::BilinearForm a(velocity, delta);
ConvectionDiffusion::LinearForm L(u0, velocity, f, delta);
                                                                                                    
// Assemble left-hand side
FEM::assemble(a, A, mesh);
                                                                                                   
// Parameters for time-stepping
real T = 0.3;
real k = 0.01;
real t = k;
                                                                                                    
// Output file
File file("temperature.pvd");

// Start time-stepping
Progress p("Time-stepping");
while ( t < T )
{
  // Assemble load vector and set boundary conditions
  FEM::assemble(L, b, mesh);
  FEM::applyBC(A, b, mesh, element, bc);
                                                                                                    
  // Solve the linear system
  solver.solve(A, x, b);
                                                                                                    
  // Save the solution to file
  file << u1;
                                                                                                    
  // Update progress
  p = t / T;
                                                                                                    
  // Move to next interval
  t += k;
  u0 = u1;
}

Note in particular the following things that you couldn't do before:

1. Function u0 = 0.0;

This initializes a Function that is constant equal to zero.

2. u0 = u1;

This updates the values of u0 to those of u1. u0 will change from
being a constant Function to a discrete Function and automatically
create a new Vector (previously x0) and copy the values from x1.
The next time, it will keep the previously created Vector and just
copy the values.

I haven't actually checked the solution so I'll just keep my fingers
crossed I didn't break anything... :-)

/Anders