← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to dolfin/fem/DirichletBC.cpp

  • Committer: Kent-Andre Mardal
  • Date: 2008-05-19 14:21:52 UTC
  • mfrom: (2668.5.1 trunk)
  • mto: (2668.1.16 trunk)
  • mto: This revision was merged to the branch mainline in revision 2670.
  • Revision ID: kent-and@simula.no-20080519142152-7zb7r4htl7111izh
merge

Show diffs side-by-side

added added

removed removed

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