← Back to branch summary

~corrado-maurini/dolfin/tao

« back to all changes in this revision

Viewing changes to dolfin/adaptivity/ErrorControl.cpp

  • Committer: corrado maurini
  • Date: 2012-12-18 12:16:08 UTC
  • mfrom: (6685.78.207 trunk)
  • Revision ID: corrado.maurini@upmc.fr-20121218121608-nk82ly9jgsld9u84
updating with trunk, fix uint in TAO solver and hacking the check for tao FindTAO.cmake

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/adaptivity/ErrorControl.cpp
75
75
  _is_linear = is_linear;
76
76
 
77
77
  // Extract and store additional function spaces
78
 
  const uint improved_dual = _residual->num_coefficients() - 1;
 
78
  const std::size_t improved_dual = _residual->num_coefficients() - 1;
79
79
  const Function& e_tmp = dynamic_cast<const Function&>(*_residual->coefficient(improved_dual));
80
80
  _E = e_tmp.function_space();
81
81
 
106
106
 
107
107
  // Extract number of coefficients in residual
108
108
  dolfin_assert(_residual);
109
 
  const uint num_coeffs = _residual->num_coefficients();
 
109
  const std::size_t num_coeffs = _residual->num_coefficients();
110
110
 
111
111
  // Attach improved dual approximation to residual
112
112
  _residual->set_coefficient(num_coeffs - 1, _Ez_h);
120
120
  }
121
121
 
122
122
  // Assemble error estimate
123
 
  info("Assembling error estimate.");
 
123
  log(PROGRESS, "Assembling error estimate.");
124
124
  const double error_estimate = assemble(*_residual);
125
125
 
126
126
  // Return estimate
130
130
void ErrorControl::compute_dual(Function& z,
131
131
   const std::vector<boost::shared_ptr<const BoundaryCondition> > bcs)
132
132
{
133
 
  info("Solving dual problem.");
 
133
  log(PROGRESS, "Solving dual problem.");
134
134
 
135
135
  // Create dual boundary conditions by homogenizing
136
136
  std::vector<boost::shared_ptr<const BoundaryCondition> > dual_bcs;
137
 
  for (uint i = 0; i < bcs.size(); i++)
 
137
  for (std::size_t i = 0; i < bcs.size(); i++)
138
138
  {
139
139
    // Only handle DirichletBCs
140
140
    const DirichletBC* bc_ptr = dynamic_cast<const DirichletBC*>(bcs[i].get());
166
166
void ErrorControl::compute_extrapolation(const Function& z,
167
167
   const std::vector<boost::shared_ptr<const BoundaryCondition> > bcs)
168
168
{
169
 
  info("Extrapolating dual solution.");
 
169
  log(PROGRESS, "Extrapolating dual solution.");
170
170
 
171
171
  // Extrapolate
172
172
  dolfin_assert(_E);
187
187
  // Create SpecialFacetFunction for the strong facet residual (R_dT)
188
188
  dolfin_assert(_a_R_dT);
189
189
  std::vector<Function> f_e;
190
 
  for (uint i = 0; i <= _R_T->geometric_dimension(); i++)
 
190
  for (std::size_t i = 0; i <= _R_T->geometric_dimension(); i++)
191
191
    f_e.push_back(Function(_a_R_dT->function_space(1)));
192
192
 
193
193
  if (f_e[0].value_rank() == 0)
234
234
  // Convert DG_0 vector to mesh function over cells
235
235
  for (CellIterator cell(mesh); !cell.end(); ++cell)
236
236
  {
237
 
    const std::vector<uint>& dofs = dofmap.cell_dofs(cell->index());
 
237
    const std::vector<DolfinIndex>& dofs = dofmap.cell_dofs(cell->index());
238
238
    dolfin_assert(dofs.size() == 1);
239
239
    indicators[cell->index()] = x[dofs[0]];
240
240
  }
245
245
                                           SpecialFacetFunction& R_dT,
246
246
                                           const Function& u)
247
247
{
248
 
  begin("Computing residual representation.");
 
248
  begin(PROGRESS, "Computing residual representation.");
249
249
 
250
250
  // Compute cell residual
251
251
  Timer timer("Computation of residual representation");
260
260
//-----------------------------------------------------------------------------
261
261
void ErrorControl::compute_cell_residual(Function& R_T, const Function& u)
262
262
{
263
 
  begin("Computing cell residual representation.");
 
263
  begin(PROGRESS, "Computing cell residual representation.");
264
264
 
265
265
  dolfin_assert(_a_R_T);
266
266
  dolfin_assert(_L_R_T);
267
267
  dolfin_assert(_cell_bubble);
268
268
 
269
269
  // Attach cell bubble to _a_R_T and _L_R_T
270
 
  const uint num_coeffs = _L_R_T->num_coefficients();
 
270
  const std::size_t num_coeffs = _L_R_T->num_coefficients();
271
271
  _a_R_T->set_coefficient(0, _cell_bubble);
272
272
  _L_R_T->set_coefficient(num_coeffs - 1, _cell_bubble);
273
273
 
292
292
 
293
293
  // Define matrices for cell-residual problems
294
294
  dolfin_assert(V.element());
295
 
  const uint N = V.element()->space_dimension();
 
295
  const std::size_t N = V.element()->space_dimension();
296
296
  arma::mat A(N, N);
297
297
  arma::mat b(N, 1);
298
298
  arma::vec x(N);
299
299
 
300
300
  // Extract cell_domains etc from right-hand side form
301
 
  const MeshFunction<uint>*
 
301
  const MeshFunction<std::size_t>*
302
302
    cell_domains = _L_R_T->cell_domains_shared_ptr().get();
303
 
  const MeshFunction<uint>*
 
303
  const MeshFunction<std::size_t>*
304
304
    exterior_facet_domains = _L_R_T->exterior_facet_domains_shared_ptr().get();
305
 
  const MeshFunction<uint>*
 
305
  const MeshFunction<std::size_t>*
306
306
    interior_facet_domains = _L_R_T->interior_facet_domains_shared_ptr().get();
307
307
 
308
308
  // Assemble and solve local linear systems
318
318
    x = arma::solve(A, b);
319
319
 
320
320
    // Get local-to-global dof map for cell
321
 
    const std::vector<uint>& dofs = dofmap.cell_dofs(cell->index());
 
321
    const std::vector<DolfinIndex>& dofs = dofmap.cell_dofs(cell->index());
322
322
 
323
323
    // Plug local solution into global vector
324
324
    dolfin_assert(R_T.vector());
331
331
                                          const Function& u,
332
332
                                          const Function& R_T)
333
333
{
334
 
  begin("Computing facet residual representation.");
 
334
  begin(PROGRESS, "Computing facet residual representation.");
335
335
 
336
336
  // Extract function space for facet residual approximation
337
337
  dolfin_assert(R_dT[0].function_space());
338
338
  const FunctionSpace& V = *R_dT[0].function_space();
339
339
  dolfin_assert(V.element());
340
 
  const uint N = V.element()->space_dimension();
 
340
  const std::size_t N = V.element()->space_dimension();
341
341
 
342
342
  // Extract mesh
343
343
  dolfin_assert(V.mesh());
351
351
  // Extract number of coefficients on right-hand side (for use with
352
352
  // attaching coefficients)
353
353
  dolfin_assert(_L_R_dT);
354
 
  const uint L_R_dT_num_coefficients = _L_R_dT->num_coefficients();
 
354
  const std::size_t L_R_dT_num_coefficients = _L_R_dT->num_coefficients();
355
355
 
356
356
  // Attach primal approximation if linear (already attached
357
357
  // otherwise).
375
375
  arma::vec x(N);
376
376
 
377
377
  // Variables to be used for the construction of the cone function
378
 
  const uint num_cells = mesh.num_cells();
 
378
  const std::size_t num_cells = mesh.num_cells();
379
379
  const std::vector<double> ones(num_cells, 1.0);
380
 
  std::vector<uint> facet_dofs(num_cells);
 
380
  std::vector<DolfinIndex> facet_dofs(num_cells);
381
381
 
382
382
  // Extract cell_domains etc from right-hand side form
383
383
  dolfin_assert(_L_R_T);
384
 
  const MeshFunction<uint>*
 
384
  const MeshFunction<std::size_t>*
385
385
    cell_domains = _L_R_T->cell_domains_shared_ptr().get();
386
 
  const MeshFunction<uint>*
 
386
  const MeshFunction<std::size_t>*
387
387
    exterior_facet_domains = _L_R_T->exterior_facet_domains_shared_ptr().get();
388
 
  const MeshFunction<uint>*
 
388
  const MeshFunction<std::size_t>*
389
389
    interior_facet_domains = _L_R_T->interior_facet_domains_shared_ptr().get();
390
390
 
391
391
  dolfin_assert(_a_R_dT);
397
397
    dolfin_assert(_cell_cone->vector());
398
398
    *(_cell_cone->vector()) = 0.0;
399
399
    facet_dofs.clear();
400
 
    const uint local_facet_dof = local_cone_dim - (dim + 1) + local_facet;
 
400
    const std::size_t local_facet_dof = local_cone_dim - (dim + 1) + local_facet;
401
401
    dolfin_assert(_cell_cone->function_space());
402
402
    dolfin_assert(_cell_cone->function_space()->dofmap());
403
403
    const GenericDofMap& cone_dofmap(*(_cell_cone->function_space()->dofmap()));
404
 
    for (uint k = 0; k < num_cells; k++)
 
404
    for (std::size_t k = 0; k < num_cells; k++)
405
405
      facet_dofs.push_back(cone_dofmap.cell_dofs(k)[local_facet_dof]);
406
406
    _cell_cone->vector()->set(&ones[0], num_cells, &facet_dofs[0]);
407
407
 
423
423
                               exterior_facet_domains, interior_facet_domains);
424
424
 
425
425
      // Non-singularize local matrix
426
 
      for (uint i = 0; i < N; ++i)
 
426
      for (std::size_t i = 0; i < N; ++i)
427
427
      {
428
428
        if (std::abs(A(i, i)) < 1.0e-10)
429
429
        {
436
436
      x = arma::solve(A, b);
437
437
 
438
438
      // Get local-to-global dof map for cell
439
 
      const std::vector<uint>& dofs = dofmap.cell_dofs(cell->index());
 
439
      const std::vector<DolfinIndex>& dofs = dofmap.cell_dofs(cell->index());
440
440
 
441
441
      // Plug local solution into global vector
442
442
      dolfin_assert(R_dT[local_facet].vector());
451
451
{
452
452
  // Create boundary conditions for extrapolated dual, and apply
453
453
  // these.
454
 
  for (uint i = 0; i < bcs.size(); i++)
 
454
  for (std::size_t i = 0; i < bcs.size(); i++)
455
455
  {
456
456
    // Only handle DirichletBCs
457
457
    const DirichletBC* bc = dynamic_cast<const DirichletBC*>(bcs[i].get());
459
459
 
460
460
    // Extract SubSpace component
461
461
    dolfin_assert(bc->function_space());
462
 
    const std::vector<uint> component = bc->function_space()->component();
 
462
    const std::vector<std::size_t> component = bc->function_space()->component();
463
463
 
464
464
    // Extract sub-domain
465
465
    boost::shared_ptr<const SubDomain> sub_domain = bc->user_sub_domain();