~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to dolfin/fem/DirichletBC.cpp

  • Committer: Anders Logg
  • Date: 2008-05-21 22:42:48 UTC
  • mfrom: (2668.6.1 trunk)
  • mto: (2668.8.3 trunk)
  • mto: This revision was merged to the branch mainline in revision 2670.
  • Revision ID: logg@simula.no-20080521224248-7baydkw3uy323fur
merge

Show diffs side-by-side

added added

removed removed

Lines of Context:
4
4
// Modified by Kristian Oelgaard, 2007
5
5
//
6
6
// First added:  2007-04-10
7
 
// Last changed: 2008-05-18
 
7
// Last changed: 2008-05-21
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>
12
13
#include <dolfin/mesh/Vertex.h>
13
14
#include <dolfin/mesh/Cell.h>
14
15
#include <dolfin/mesh/Facet.h>
32
33
  : BoundaryCondition(), g(g), _mesh(mesh),
33
34
    sub_domains(0), sub_domain(0), sub_domains_local(false),
34
35
    method(method), user_sub_domain(&sub_domain)
35
 
 
36
36
{
37
37
  // Initialize sub domain markers
38
 
  init(sub_domain);
 
38
  initFromSubDomain(sub_domain);
39
39
}
40
40
//-----------------------------------------------------------------------------
41
41
DirichletBC::DirichletBC(Function& g,
51
51
//-----------------------------------------------------------------------------
52
52
DirichletBC::DirichletBC(Function& g,
53
53
                         Mesh& mesh,
 
54
                         uint sub_domain,
 
55
                         BCMethod method)
 
56
  : BoundaryCondition(), g(g), _mesh(mesh),
 
57
    sub_domains(0), sub_domain(sub_domain), sub_domains_local(false),
 
58
    method(method), user_sub_domain(0)
 
59
{
 
60
  // Initialize sub domain markers
 
61
  initFromMesh();
 
62
}
 
63
//-----------------------------------------------------------------------------
 
64
DirichletBC::DirichletBC(Function& g,
 
65
                         Mesh& mesh,
54
66
                         SubDomain& sub_domain,
55
67
                         const SubSystem& sub_system,
56
68
                         BCMethod method)
58
70
    sub_domains(0), sub_domain(0), sub_domains_local(false),
59
71
    sub_system(sub_system), method(method), user_sub_domain(&sub_domain)
60
72
{
61
 
  // Set sub domain markers
62
 
  init(sub_domain);
 
73
  // Initialize sub domain markers
 
74
  initFromSubDomain(sub_domain);
63
75
}
64
76
//-----------------------------------------------------------------------------
65
77
DirichletBC::DirichletBC(Function& g,
76
88
//-----------------------------------------------------------------------------
77
89
DirichletBC::DirichletBC(Function& g,
78
90
                         Mesh& mesh,
 
91
                         uint sub_domain,
 
92
                         const SubSystem& sub_system,
79
93
                         BCMethod method)
80
94
  : BoundaryCondition(), g(g), _mesh(mesh),
81
 
    sub_domains(0), sub_domain(0), sub_domains_local(false),
82
 
    method(method), user_sub_domain(0)
 
95
    sub_domains(), sub_domain(sub_domain), sub_domains_local(false),
 
96
    sub_system(sub_system), method(method), user_sub_domain(0)
83
97
{
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
98
  // Initialize sub domain markers
97
 
  init(sub_domain);
 
99
  initFromMesh();
98
100
}
99
101
//-----------------------------------------------------------------------------
100
102
DirichletBC::~DirichletBC()
174
176
    Progress p(s, _mesh.size(D - 1));
175
177
    for (FacetIterator facet(_mesh); !facet.end(); ++facet)
176
178
    {
177
 
      // Skip facets not inside the sub domain
 
179
      // Skip facets not inside the sub domain      
178
180
      if ((*sub_domains)(*facet) != sub_domain)
179
181
      {
180
182
        p++;
236
238
  return _mesh;
237
239
}
238
240
//-----------------------------------------------------------------------------
239
 
void DirichletBC::init(SubDomain& sub_domain)
 
241
void DirichletBC::initFromSubDomain(SubDomain& sub_domain)
240
242
{
241
243
  cout << "Creating sub domain markers for boundary condition." << endl;
242
244
 
252
254
  sub_domain.mark(*sub_domains, 0);
253
255
}
254
256
//-----------------------------------------------------------------------------
 
257
void DirichletBC::initFromMesh()
 
258
{
 
259
  cout << "Creating sub domain markers for boundary condition." << endl;
 
260
  
 
261
  // Get data
 
262
  Array<uint>* facet_cells   = _mesh.data().array("boundary facet cells");
 
263
  Array<uint>* facet_numbers = _mesh.data().array("boundary facet numbers");
 
264
  Array<uint>* indicators    = _mesh.data().array("boundary indicators");
 
265
 
 
266
  // Check data
 
267
  if (!facet_cells)
 
268
    error("Mesh data \"boundary facet cells\" not available.");
 
269
  if (!facet_numbers)
 
270
    error("Mesh data \"boundary facet numbers\" not available.");
 
271
  if (!indicators)
 
272
    error("Mesh data \"boundary indicators\" not available.");
 
273
 
 
274
  // Get size
 
275
  const uint size = facet_cells->size();
 
276
  dolfin_assert(size == facet_numbers->size());
 
277
  dolfin_assert(size == indicators->size());
 
278
 
 
279
  // Create mesh function for sub domain markers on facets
 
280
  const uint dim = _mesh.topology().dim();
 
281
  _mesh.init(dim - 1);
 
282
  sub_domains = new MeshFunction<uint>(_mesh, dim - 1);
 
283
  sub_domains_local = true;
 
284
 
 
285
  // Compute the maximum boundary indicator
 
286
  uint maxid = 0;
 
287
  for (uint i = 0; i < size; i++)
 
288
    maxid = std::max(maxid, (*indicators)[i]);
 
289
  cout << "maxid = " << maxid << endl;
 
290
 
 
291
  // Mark everything with the maximum value + 1
 
292
  (*sub_domains) = maxid + 1;
 
293
 
 
294
  // Mark facets according to data
 
295
  for (uint i = 0; i < size; i++)
 
296
  {
 
297
    // Get cell incident with facet
 
298
    Cell cell(_mesh, (*facet_cells)[i]);
 
299
 
 
300
    // Get facet
 
301
    uint facet = cell.entities(dim - 1)[(*facet_numbers)[i]];
 
302
 
 
303
    // Set marker
 
304
    sub_domains->set(facet, (*indicators)[i]);
 
305
  }
 
306
}
 
307
//-----------------------------------------------------------------------------
255
308
void DirichletBC::computeBCTopological(std::map<uint, real>& boundary_values,
256
309
                                       Facet& facet,
257
310
                                       BoundaryCondition::LocalData& data)
470
523
  std::map<uint, real>::const_iterator boundary_value;
471
524
  uint i = 0;
472
525
  for (boundary_value = boundary_values.begin(); boundary_value != boundary_values.end(); ++boundary_value)
473
 
  {
474
 
    dofs[i++]     = boundary_value->first;
475
 
  }
 
526
    dofs[i++] = boundary_value->first;
476
527
 
477
528
  // Modify linear system (A_ii = 1)
478
529
  A.zero(boundary_values.size(), dofs);