← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to dolfin/fem/DirichletBC.cpp

  • Committer: Anders Logg
  • Date: 2008-05-22 08:19:48 UTC
  • mto: (2668.9.6 trunk)
  • mto: This revision was merged to the branch mainline in revision 2670.
  • Revision ID: logg@simula.no-20080522081948-u0ur5rmsfg76yfvs
Major cleanup of DirichletBC:
 - Store boundary facets in array, common for all different representations
 - BCs may now be specified in 3 differenet ways: SubDomain, MeshFunction, or as Mesh data
 - Should speed up repeated application of BCs
 - Tested on a few different applications, but please check that nothing broke.
   In particular, geometric and pointwise application has not been tested.

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/fem/DirichletBC.cpp
31
31
                         SubDomain& sub_domain,
32
32
                         BCMethod method)
33
33
  : BoundaryCondition(), g(g), _mesh(mesh),
34
 
    sub_domains(0), sub_domain(0), sub_domains_local(false),
35
34
    method(method), user_sub_domain(&sub_domain)
36
35
{
37
 
  // Initialize sub domain markers
38
36
  initFromSubDomain(sub_domain);
39
37
}
40
38
//-----------------------------------------------------------------------------
43
41
                         uint sub_domain,
44
42
                         BCMethod method)
45
43
  : BoundaryCondition(), g(g), _mesh(sub_domains.mesh()),
46
 
    sub_domains(&sub_domains), sub_domain(sub_domain), sub_domains_local(false),
47
44
    method(method), user_sub_domain(0)
48
45
{
49
 
  // Do nothing
 
46
  initFromMeshFunction(sub_domains, sub_domain);
50
47
}
51
48
//-----------------------------------------------------------------------------
52
49
DirichletBC::DirichletBC(Function& g,
54
51
                         uint sub_domain,
55
52
                         BCMethod method)
56
53
  : BoundaryCondition(), g(g), _mesh(mesh),
57
 
    sub_domains(0), sub_domain(sub_domain), sub_domains_local(false),
58
54
    method(method), user_sub_domain(0)
59
55
{
60
 
  // Initialize sub domain markers
61
56
  initFromMesh(sub_domain);
62
57
}
63
58
//-----------------------------------------------------------------------------
67
62
                         const SubSystem& sub_system,
68
63
                         BCMethod method)
69
64
  : BoundaryCondition(), g(g), _mesh(mesh),
70
 
    sub_domains(0), sub_domain(0), sub_domains_local(false),
71
 
    sub_system(sub_system), method(method), user_sub_domain(&sub_domain)
 
65
    method(method), user_sub_domain(&sub_domain),
 
66
    sub_system(sub_system)
72
67
{
73
 
  // Initialize sub domain markers
74
68
  initFromSubDomain(sub_domain);
75
69
}
76
70
//-----------------------------------------------------------------------------
80
74
                         const SubSystem& sub_system,
81
75
                         BCMethod method)
82
76
  : BoundaryCondition(), g(g), _mesh(sub_domains.mesh()),
83
 
    sub_domains(&sub_domains), sub_domain(sub_domain), sub_domains_local(false),
84
 
    sub_system(sub_system), method(method), user_sub_domain(0)
 
77
    method(method), user_sub_domain(0),
 
78
    sub_system(sub_system)
85
79
{
86
 
  // Do nothing
 
80
  initFromMeshFunction(sub_domains, sub_domain);
87
81
}
88
82
//-----------------------------------------------------------------------------
89
83
DirichletBC::DirichletBC(Function& g,
92
86
                         const SubSystem& sub_system,
93
87
                         BCMethod method)
94
88
  : BoundaryCondition(), g(g), _mesh(mesh),
95
 
    sub_domains(), sub_domain(sub_domain), sub_domains_local(false),
96
 
    sub_system(sub_system), method(method), user_sub_domain(0)
 
89
    method(method), user_sub_domain(0),
 
90
    sub_system(sub_system)
97
91
{
98
 
  // Initialize sub domain markers
99
92
  initFromMesh(sub_domain);
100
93
}
101
94
//-----------------------------------------------------------------------------
102
95
DirichletBC::~DirichletBC()
103
96
{
104
 
  // Delete sub domain markers if created locally
105
 
  if ( sub_domains_local )
106
 
    delete sub_domains;
 
97
  // Do nothing
107
98
}
108
99
//-----------------------------------------------------------------------------
109
100
void DirichletBC::apply(GenericMatrix& A, GenericVector& b, const Form& form)
132
123
void DirichletBC::apply(GenericMatrix& A, GenericVector& b,
133
124
                        const GenericVector* x, const DofMap& dof_map, const ufc::form& form)
134
125
{
135
 
  // FIXME: How do we reuse the dof map for u?
136
 
 
137
126
  // Simple check
138
127
  const uint N = dof_map.global_dimension();
139
128
  if (N != A.size(0) || N != A.size(1))
141
130
  if (N != b.size())
142
131
    error("Incorrect dimension of matrix for application of boundary conditions. Did you assemble it on a different mesh?");
143
132
 
144
 
  // Set message string
145
 
  std::string s;
146
 
  if (method == topological)
147
 
    s = "Computing Dirichlet boundary values";
148
 
  else if (method == geometrical)
149
 
    s = "Computing Dirichlet boundary values (geometrical approach)";
150
 
  else
151
 
    s = "Computing Dirichlet boundary values (pointwise approach)";
152
 
  
153
 
  // Make sure we have the facet - cell connectivity
154
 
  const uint D = _mesh.topology().dim();
155
 
  if (method == topological)
156
 
    _mesh.init(D - 1, D);
 
133
  // A map to hold the mapping from boundary dofs to boundary values
 
134
  std::map<uint, real> boundary_values;
157
135
 
158
136
  // Create local data for application of boundary conditions
159
137
  BoundaryCondition::LocalData data(form, _mesh, dof_map, sub_system);
160
 
  
161
 
  // A map to hold the mapping from boundary dofs to boundary values
162
 
  std::map<uint, real> boundary_values;
163
 
 
164
 
  // Choose strategy
165
 
  if (method == pointwise)
166
 
  {
167
 
    Progress p(s, _mesh.size(D));
168
 
    for (CellIterator cell(_mesh); !cell.end(); ++cell)
169
 
    {
170
 
      computeBCPointwise(boundary_values, *cell, data);
171
 
      p++;
172
 
    }
173
 
  }
174
 
  else
175
 
  {
176
 
    // Iterate over the facets of the mesh
177
 
    Progress p(s, _mesh.size(D - 1));
178
 
    for (FacetIterator facet(_mesh); !facet.end(); ++facet)
179
 
    {
180
 
      // Skip facets not inside the sub domain      
181
 
      if ((*sub_domains)(*facet) != sub_domain)
182
 
      {
183
 
        p++;
184
 
        continue;
185
 
      }
186
 
 
187
 
      // Chose strategy
188
 
      if (method == topological)
189
 
        computeBCTopological(boundary_values, *facet, data);
190
 
      else
191
 
        computeBCGeometrical(boundary_values, *facet, data);
192
 
    
193
 
      // Update process
194
 
      p++;
195
 
    }
196
 
  }
 
138
 
 
139
  // Compute dofs and values
 
140
  computeBC(boundary_values, data);
197
141
 
198
142
  // Copy boundary value data to arrays
199
143
  uint* dofs = new uint[boundary_values.size()];
205
149
    dofs[i]     = boundary_value->first;
206
150
    values[i++] = boundary_value->second;
207
151
  }
208
 
 
 
152
  
209
153
  // Modify boundary values for nonlinear problems
210
154
  if (x)
211
155
  {
214
158
    for (uint i = 0; i < boundary_values.size(); i++)
215
159
      values[i] -= x_values[i];
216
160
  }
217
 
 
 
161
  
218
162
  message("Applying boundary conditions to linear system.");
219
163
  
220
164
  // Modify RHS vector (b[i] = value)
221
165
  b.set(values, boundary_values.size(), dofs);
222
 
 
 
166
  
223
167
  // Modify linear system (A_ii = 1)
224
168
  A.ident(boundary_values.size(), dofs);
225
 
 
 
169
  
226
170
  // Clear temporary arrays
227
171
  delete [] dofs;
228
172
  delete [] values;
229
 
 
 
173
  
230
174
  // Finalise changes to A
231
175
  A.apply();
232
 
 
 
176
  
233
177
  // Finalise changes to b
234
178
  b.apply();
235
179
}
236
180
//-----------------------------------------------------------------------------
 
181
void DirichletBC::zero(GenericMatrix& A, const Form& form)
 
182
{
 
183
  zero(A, form.dofMaps()[1], form.form());
 
184
}
 
185
//-----------------------------------------------------------------------------
 
186
void DirichletBC::zero(GenericMatrix& A, const DofMap& dof_map, const ufc::form& form)
 
187
{
 
188
  // Simple check
 
189
  const uint N = dof_map.global_dimension();
 
190
  if (N != A.size(0) || N != A.size(1))
 
191
    error("Incorrect dimension of matrix for application of boundary conditions. Did you assemble it on a different mesh?");
 
192
 
 
193
  // A map to hold the mapping from boundary dofs to boundary values
 
194
  std::map<uint, real> boundary_values;
 
195
 
 
196
  // Create local data for application of boundary conditions
 
197
  BoundaryCondition::LocalData data(form, _mesh, dof_map, sub_system);
 
198
 
 
199
  // Compute dofs and values
 
200
  computeBC(boundary_values, data);
 
201
 
 
202
  // Copy boundary value data to arrays
 
203
  uint* dofs = new uint[boundary_values.size()];
 
204
  std::map<uint, real>::const_iterator boundary_value;
 
205
  uint i = 0;
 
206
  for (boundary_value = boundary_values.begin(); boundary_value != boundary_values.end(); ++boundary_value)
 
207
    dofs[i++] = boundary_value->first;
 
208
 
 
209
  // Modify linear system (A_ii = 1)
 
210
  A.zero(boundary_values.size(), dofs);
 
211
 
 
212
  // Finalise changes to A
 
213
  A.apply();
 
214
 
 
215
  // Clear temporary arrays
 
216
  delete [] dofs;
 
217
}
 
218
//-----------------------------------------------------------------------------
237
219
Mesh& DirichletBC::mesh()
238
220
{
239
221
  return _mesh;
241
223
//-----------------------------------------------------------------------------
242
224
void DirichletBC::initFromSubDomain(SubDomain& sub_domain)
243
225
{
244
 
  cout << "Creating sub domain markers for boundary condition." << endl;
 
226
  dolfin_assert(facets.size() == 0);
 
227
 
 
228
  // FIXME: This can be made more efficient, we should be able to
 
229
  // FIXME: extract the facets without first creating a MeshFunction on
 
230
  // FIXME: the entire mesh and then extracting the subset. This is done
 
231
  // FIXME: mainly for convenience (we may reuse mark() in SubDomain).
245
232
 
246
233
  // Create mesh function for sub domain markers on facets
247
 
  _mesh.init(_mesh.topology().dim() - 1);
248
 
  sub_domains = new MeshFunction<uint>(_mesh, _mesh.topology().dim() - 1);
249
 
  sub_domains_local = true;
 
234
  const uint dim = _mesh.topology().dim();
 
235
  _mesh.init(dim - 1);
 
236
  MeshFunction<uint> sub_domains(_mesh, dim - 1);
250
237
 
251
238
  // Mark everything as sub domain 1
252
 
  (*sub_domains) = 1;
 
239
  sub_domains = 1;
253
240
  
254
241
  // Mark the sub domain as sub domain 0
255
 
  sub_domain.mark(*sub_domains, 0);
 
242
  sub_domain.mark(sub_domains, 0);
 
243
 
 
244
  // Initialize from mesh function
 
245
  initFromMeshFunction(sub_domains, 0);
 
246
}
 
247
//-----------------------------------------------------------------------------
 
248
void DirichletBC::initFromMeshFunction(MeshFunction<uint>& sub_domains,
 
249
                                       uint sub_domain)
 
250
{
 
251
  dolfin_assert(facets.size() == 0);
 
252
 
 
253
  // Make sure we have the facet - cell connectivity
 
254
  const uint dim = _mesh.topology().dim();
 
255
  _mesh.init(dim - 1, dim);
 
256
 
 
257
  // Build set of boundary facets
 
258
  for (FacetIterator facet(_mesh); !facet.end(); ++facet)
 
259
  {
 
260
    // Skip facets not on this boundary
 
261
    if (sub_domains(*facet) != sub_domain)
 
262
      continue;
 
263
 
 
264
    // Get cell to which facet belongs (there may be two, but pick first)
 
265
    Cell cell(_mesh, facet->entities(dim)[0]);
 
266
    
 
267
    // Get local index of facet with respect to the cell
 
268
    const uint facet_number = cell.index(*facet);
 
269
 
 
270
    // Copy data
 
271
    facets.push_back(std::pair<uint, uint>(cell.index(), facet_number));
 
272
  }
256
273
}
257
274
//-----------------------------------------------------------------------------
258
275
void DirichletBC::initFromMesh(uint sub_domain)
291
308
  }
292
309
}
293
310
//-----------------------------------------------------------------------------
 
311
void DirichletBC::computeBC(std::map<uint, real>& boundary_values,
 
312
                            BoundaryCondition::LocalData& data)
 
313
{
 
314
  // Choose strategy
 
315
  switch (method)
 
316
  {
 
317
  case topological:
 
318
    computeBCTopological(boundary_values, data);
 
319
    break;
 
320
  case geometric:
 
321
    computeBCGeometric(boundary_values, data);
 
322
    break;
 
323
  case pointwise:
 
324
    computeBCPointwise(boundary_values, data);
 
325
    break;
 
326
  default:
 
327
    error("Unknown method for application of boundary conditions.");
 
328
  }
 
329
}
 
330
//-----------------------------------------------------------------------------
294
331
void DirichletBC::computeBCTopological(std::map<uint, real>& boundary_values,
295
 
                                       Facet& facet,
296
332
                                       BoundaryCondition::LocalData& data)
297
333
{
298
 
  // Get cell to which facet belongs (there may be two, but pick first)
299
 
  Cell cell(_mesh, facet.entities(facet.dim() + 1)[0]);
300
 
  UFCCell ufc_cell(cell);
301
 
  
302
 
  // Get local index of facet with respect to the cell
303
 
  const uint local_facet = cell.index(facet);
304
 
  
305
 
  // Interpolate function on cell
306
 
  g.interpolate(data.w, ufc_cell, *data.finite_element, cell, local_facet);
307
 
  
308
 
  // Tabulate dofs on cell
309
 
  data.dof_map->tabulate_dofs(data.cell_dofs, ufc_cell);
310
 
 
311
 
  // Tabulate which dofs are on the facet
312
 
  data.dof_map->tabulate_facet_dofs(data.facet_dofs, local_facet);
313
 
  
314
 
  // Debugging print:
315
 
  /* 
316
 
  cout << endl << "Handling BC's for:" << endl;
317
 
  cout << "Cell:  " << facet.entities(facet.dim() + 1)[0] << endl;
318
 
  cout << "Facet: " << local_facet << endl;
319
 
  */
320
 
  
321
 
  // Pick values for facet
322
 
  for (uint i = 0; i < data.dof_map->num_facet_dofs(); i++)
 
334
  // Iterate over facets
 
335
  Progress p("Computing Dirichlet boundary values, topological search", facets.size());
 
336
  for (uint f = 0; f < facets.size(); f++)
323
337
  {
324
 
    const uint dof = data.offset + data.cell_dofs[data.facet_dofs[i]];
325
 
    const real value = data.w[data.facet_dofs[i]];
326
 
    boundary_values[dof] = value;
327
 
    //cout << "Setting BC value: i = " << i << ", dof = " << dof << ", value = " << value << endl;
 
338
    // Get cell number and local facet number
 
339
    uint cell_number = facets[f].first;
 
340
    uint facet_number = facets[f].second;
 
341
 
 
342
    // Create cell
 
343
    Cell cell(_mesh, cell_number);
 
344
    UFCCell ufc_cell(cell);
 
345
 
 
346
    // Interpolate function on cell
 
347
    g.interpolate(data.w, ufc_cell, *data.finite_element, cell, facet_number);
 
348
    
 
349
    // Tabulate dofs on cell
 
350
    data.dof_map->tabulate_dofs(data.cell_dofs, ufc_cell);
 
351
    
 
352
    // Tabulate which dofs are on the facet
 
353
    data.dof_map->tabulate_facet_dofs(data.facet_dofs, facet_number);
 
354
    
 
355
    // Debugging print:
 
356
    /* 
 
357
       cout << endl << "Handling BC's for:" << endl;
 
358
       cout << "Cell:  " << facet.entities(facet.dim() + 1)[0] << endl;
 
359
       cout << "Facet: " << local_facet << endl;
 
360
    */
 
361
    
 
362
    // Pick values for facet
 
363
    for (uint i = 0; i < data.dof_map->num_facet_dofs(); i++)
 
364
    {
 
365
      const uint dof = data.offset + data.cell_dofs[data.facet_dofs[i]];
 
366
      const real value = data.w[data.facet_dofs[i]];
 
367
      boundary_values[dof] = value;
 
368
      //cout << "Setting BC value: i = " << i << ", dof = " << dof << ", value = " << value << endl;
 
369
    }
 
370
 
 
371
    p++;
328
372
  }
329
 
 
330
373
}
331
374
//-----------------------------------------------------------------------------
332
 
void DirichletBC::computeBCGeometrical(std::map<uint, real>& boundary_values,
333
 
                                       Facet& facet,
334
 
                                       BoundaryCondition::LocalData& data)
 
375
void DirichletBC::computeBCGeometric(std::map<uint, real>& boundary_values,
 
376
                                     BoundaryCondition::LocalData& data)
335
377
{
336
 
  // Loop the vertices associated with the facet
337
 
  for (VertexIterator vertex(facet); !vertex.end(); ++vertex)
 
378
  // Initialize facets, needed for geometric search
 
379
  message("Computing facets, needed for geometric application of boundary conditions.");
 
380
  _mesh.init(_mesh.topology().dim() - 1);
 
381
 
 
382
  // Iterate over facets
 
383
  Progress p("Computing Dirichlet boundary values, geometric search", facets.size());
 
384
  for (uint f = 0; f < facets.size(); f++)
338
385
  {
339
 
    // Loop the cells associated with the vertex
340
 
    for (CellIterator c(*vertex); !c.end(); ++c)
 
386
    // Get cell number and local facet number
 
387
    uint cell_number = facets[f].first;
 
388
    uint facet_number = facets[f].second;
 
389
 
 
390
    // Create facet
 
391
    Cell cell(_mesh, cell_number);
 
392
    Facet facet(_mesh, cell.entities(_mesh.topology().dim() - 1)[facet_number]);
 
393
 
 
394
    // Loop the vertices associated with the facet
 
395
    for (VertexIterator vertex(facet); !vertex.end(); ++vertex)
341
396
    {
342
 
      UFCCell ufc_cell(*c);
343
 
      
344
 
      // Interpolate function on cell
345
 
      g.interpolate(data.w, ufc_cell, *data.finite_element, *c);
346
 
      
347
 
      // Tabulate dofs on cell, and their coordinates
348
 
      data.dof_map->tabulate_dofs(data.cell_dofs, ufc_cell);
349
 
      data.dof_map->tabulate_coordinates(data.coordinates, ufc_cell);
350
 
      
351
 
      // Loop all dofs on cell
352
 
      for (uint i = 0; i < data.dof_map->local_dimension(); ++i)
 
397
      // Loop the cells associated with the vertex
 
398
      for (CellIterator c(*vertex); !c.end(); ++c)
353
399
      {
354
 
        // Check if the coordinates are on current facet and thus on boundary
355
 
        if (!onFacet(data.coordinates[i], facet))
356
 
          continue;
357
 
        
358
 
        // Set boundary value
359
 
        const uint dof = data.offset + data.cell_dofs[i];
360
 
        const real value = data.w[i];
361
 
        boundary_values[dof] = value;
 
400
        UFCCell ufc_cell(*c);
 
401
        
 
402
        // Interpolate function on cell
 
403
        g.interpolate(data.w, ufc_cell, *data.finite_element, *c);
 
404
        
 
405
        // Tabulate dofs on cell, and their coordinates
 
406
        data.dof_map->tabulate_dofs(data.cell_dofs, ufc_cell);
 
407
        data.dof_map->tabulate_coordinates(data.coordinates, ufc_cell);
 
408
        
 
409
        // Loop all dofs on cell
 
410
        for (uint i = 0; i < data.dof_map->local_dimension(); ++i)
 
411
        {
 
412
          // Check if the coordinates are on current facet and thus on boundary
 
413
          if (!onFacet(data.coordinates[i], facet))
 
414
            continue;
 
415
          
 
416
          // Set boundary value
 
417
          const uint dof = data.offset + data.cell_dofs[i];
 
418
          const real value = data.w[i];
 
419
          boundary_values[dof] = value;
 
420
        }
362
421
      }
363
422
    }
364
423
  }
365
424
}
366
425
//-----------------------------------------------------------------------------
367
426
void DirichletBC::computeBCPointwise(std::map<uint, real>& boundary_values,
368
 
                                     Cell& cell,
369
427
                                     BoundaryCondition::LocalData& data)
370
428
{
371
429
  dolfin_assert(user_sub_domain);
372
430
 
373
 
  // Interpolate function on cell
374
 
  UFCCell ufc_cell(cell);
375
 
  g.interpolate(data.w, ufc_cell, *data.finite_element, cell);
376
 
      
377
 
  // Tabulate dofs on cell, and their coordinates
378
 
  data.dof_map->tabulate_dofs(data.cell_dofs, ufc_cell);
379
 
  data.dof_map->tabulate_coordinates(data.coordinates, ufc_cell);
380
 
      
381
 
  // Loop all dofs on cell
382
 
  for (uint i = 0; i < data.dof_map->local_dimension(); ++i)
 
431
  // Iterate over cells
 
432
  Progress p("Computing Dirichlet boundary values, pointwise search", facets.size());
 
433
  for (CellIterator cell(_mesh); !cell.end(); ++cell)
383
434
  {
384
 
    // Check if the coordinates are part of the sub domain
385
 
    if ( !user_sub_domain->inside(data.coordinates[i], false) )
386
 
      continue;
 
435
    // Interpolate function on cell
 
436
    UFCCell ufc_cell(*cell);
 
437
    g.interpolate(data.w, ufc_cell, *data.finite_element, *cell);
 
438
    
 
439
    // Tabulate dofs on cell, and their coordinates
 
440
    data.dof_map->tabulate_dofs(data.cell_dofs, ufc_cell);
 
441
    data.dof_map->tabulate_coordinates(data.coordinates, ufc_cell);
 
442
    
 
443
    // Loop all dofs on cell
 
444
    for (uint i = 0; i < data.dof_map->local_dimension(); ++i)
 
445
    {
 
446
      // Check if the coordinates are part of the sub domain
 
447
      if ( !user_sub_domain->inside(data.coordinates[i], false) )
 
448
        continue;
 
449
      
 
450
      // Set boundary value
 
451
      const uint dof = data.offset + data.cell_dofs[i];
 
452
      const real value = data.w[i];
 
453
      boundary_values[dof] = value;
 
454
    }
387
455
 
388
 
    // Set boundary value
389
 
    const uint dof = data.offset + data.cell_dofs[i];
390
 
    const real value = data.w[i];
391
 
    boundary_values[dof] = value;
 
456
    p++;
392
457
  }
393
458
}
394
459
//-----------------------------------------------------------------------------
444
509
  return false;
445
510
}
446
511
//-----------------------------------------------------------------------------
447
 
void DirichletBC::zero(GenericMatrix& A, const Form& form)
448
 
{
449
 
  zero(A, form.dofMaps()[1], form.form());
450
 
}
451
 
//-----------------------------------------------------------------------------
452
 
void DirichletBC::zero(GenericMatrix& A, const DofMap& dof_map, const ufc::form& form)
453
 
{
454
 
  // FIXME: How do we reuse the dof map for u?
455
 
 
456
 
  std::string s;
457
 
  if (method == topological)
458
 
    s = "Applying Dirichlet boundary conditions to linear system";
459
 
  else if (method == geometrical)
460
 
    s = "Applying Dirichlet boundary conditions to linear system (geometrical approach)";
461
 
  else
462
 
    s = "Applying Dirichlet boundary conditions to linear system (pointwise approach)";
463
 
  
464
 
  // Make sure we have the facet - cell connectivity
465
 
  const uint D = _mesh.topology().dim();
466
 
  if (method == topological)
467
 
    _mesh.init(D - 1, D);
468
 
 
469
 
  // Create local data for application of boundary conditions
470
 
  BoundaryCondition::LocalData data(form, _mesh, dof_map, sub_system);
471
 
  
472
 
  // A map to hold the mapping from boundary dofs to boundary values
473
 
  std::map<uint, real> boundary_values;
474
 
 
475
 
  if (method == pointwise)
476
 
  {
477
 
    Progress p(s, _mesh.size(D));
478
 
    for (CellIterator cell(_mesh); !cell.end(); ++cell)
479
 
    {
480
 
      computeBCPointwise(boundary_values, *cell, data);
481
 
      p++;
482
 
    }
483
 
  }
484
 
  else
485
 
  {
486
 
    // Iterate over the facets of the mesh
487
 
    Progress p("Applying Dirichlet boundary conditions", _mesh.size(D - 1));
488
 
    for (FacetIterator facet(_mesh); !facet.end(); ++facet)
489
 
    {
490
 
      // Skip facets not inside the sub domain
491
 
      if ((*sub_domains)(*facet) != sub_domain)
492
 
      {
493
 
        p++;
494
 
        continue;
495
 
      }
496
 
 
497
 
      // Chose strategy
498
 
      if (method == topological)
499
 
        computeBCTopological(boundary_values, *facet, data);
500
 
      else
501
 
        computeBCGeometrical(boundary_values, *facet, data);
502
 
    
503
 
      // Update process
504
 
      p++;
505
 
    }
506
 
  }
507
 
 
508
 
  // Copy boundary value data to arrays
509
 
  uint* dofs = new uint[boundary_values.size()];
510
 
  std::map<uint, real>::const_iterator boundary_value;
511
 
  uint i = 0;
512
 
  for (boundary_value = boundary_values.begin(); boundary_value != boundary_values.end(); ++boundary_value)
513
 
    dofs[i++] = boundary_value->first;
514
 
 
515
 
  // Modify linear system (A_ii = 1)
516
 
  A.zero(boundary_values.size(), dofs);
517
 
 
518
 
  // Finalise changes to A
519
 
  A.apply();
520
 
 
521
 
  // Clear temporary arrays
522
 
  delete [] dofs;
523
 
}
524
 
//-----------------------------------------------------------------------------
525
512
void DirichletBC::setSubSystem(SubSystem sub_system)
526
513
{
527
514
  this->sub_system = sub_system;