← Back to branch summary

~ubuntu-branches/ubuntu/saucy/deal.ii/saucy

« back to all changes in this revision

Viewing changes to deal.II/source/multigrid/mg_transfer_prebuilt.cc

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV, Adam C. Powell, IV, Denis Barbier
  • Date: 2010-07-29 13:47:01 UTC
  • mfrom: (3.1.3 sid)
  • Revision ID: james.westby@ubuntu.com-20100729134701-akb8jb3stwge8tcm
Tags: 6.3.1-1
http://bugs.debian.org/562332
[ Adam C. Powell, IV ]
* Changed to source format 3.0 (quilt).
* Changed maintainer to debian-science with Adam Powell as uploader.
* Added source lintian overrides about Adam Powell's name.
* Added Vcs info on git repository.
* Bumped Standards-Version.
* Changed stamp-patch to patch target and fixed its application criterion.
* Moved make_dependencies and expand_instantiations to a versioned directory
  to avoid shlib package conflicts.

[ Denis Barbier ]
* New upstream release (closes: #562332).
  + Added libtbb support.
  + Forward-ported all patches.
* Updates for new PETSc version, including workaround for different versions
  of petsc and slepc.
* Add debian/watch.
* Update to debhelper 7.
* Added pdebuild patch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
deal.II/source/multigrid/mg_transfer_prebuilt.cc
1
1
//---------------------------------------------------------------------------
2
 
//    $Id: mg_transfer_prebuilt.cc 18724 2009-04-23 23:06:37Z bangerth $
 
2
//    $Id: mg_transfer_prebuilt.cc 20949 2010-04-06 14:42:10Z bangerth $
3
3
//    Version: $Name$
4
4
//
5
 
//    Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
 
5
//    Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
6
6
//
7
7
//    This file is subject to QPL and may not be  distributed
8
8
//    without copyright and license information. Please refer
12
12
//---------------------------------------------------------------------------
13
13
 
14
14
#include <base/logstream.h>
 
15
#include <base/function.h>
15
16
 
16
17
#include <lac/vector.h>
17
18
#include <lac/block_vector.h>
23
24
#include <fe/fe.h>
24
25
#include <multigrid/mg_dof_handler.h>
25
26
#include <multigrid/mg_dof_accessor.h>
 
27
#include <multigrid/mg_tools.h>
26
28
#include <multigrid/mg_transfer.h>
27
29
#include <multigrid/mg_transfer.templates.h>
28
30
 
29
31
DEAL_II_NAMESPACE_OPEN
30
32
 
31
 
 
32
 
 
33
 
template <typename number>
 
33
namespace
 
34
{
 
35
                                     /**
 
36
                                      * Adjust vectors on all levels to
 
37
                                      * correct size.  Here, we just
 
38
                                      * count the numbers of degrees
 
39
                                      * of freedom on each level and
 
40
                                      * @p reinit each level vector
 
41
                                      * to this length.
 
42
                                      * For compatibility reasons with
 
43
                                      * the next function
 
44
                                      * the target_component is added
 
45
                                      * here but is not used.
 
46
                                      */
 
47
  template <int dim, typename number, int spacedim>
 
48
  void
 
49
  reinit_vector (const dealii::MGDoFHandler<dim,spacedim> &mg_dof,
 
50
                 std::vector<unsigned int> ,
 
51
                 MGLevelObject<dealii::Vector<number> > &v)
 
52
  {
 
53
    for (unsigned int level=v.get_minlevel();
 
54
         level<=v.get_maxlevel();++level)
 
55
      {
 
56
        unsigned int n = mg_dof.n_dofs (level);
 
57
        v[level].reinit(n);
 
58
      }
 
59
 
 
60
  }
 
61
 
 
62
 
 
63
                                     /**
 
64
                                      * Adjust vectors on all levels to
 
65
                                      * correct size.  Here, we just
 
66
                                      * count the numbers of degrees
 
67
                                      * of freedom on each level and
 
68
                                      * @p reinit each level vector
 
69
                                      * to this length. The target_component
 
70
                                      * is handed to MGTools::count_dofs_per_block.
 
71
                                      * See for documentation there.
 
72
                                      */
 
73
  template <int dim, typename number, int spacedim>
 
74
  void
 
75
  reinit_vector (const dealii::MGDoFHandler<dim,spacedim> &mg_dof,
 
76
                 std::vector<unsigned int> target_component,
 
77
                 MGLevelObject<BlockVector<number> > &v)
 
78
  {
 
79
    const unsigned int n_blocks = mg_dof.get_fe().n_blocks();
 
80
    if (target_component.size()==0)
 
81
      {
 
82
        target_component.resize(n_blocks);
 
83
        for (unsigned int i=0;i<n_blocks;++i)
 
84
          target_component[i] = i;
 
85
      }
 
86
    Assert(target_component.size()==n_blocks,
 
87
           ExcDimensionMismatch(target_component.size(),n_blocks));
 
88
    const unsigned int max_block
 
89
      = *std::max_element (target_component.begin(),
 
90
                           target_component.end());
 
91
    const unsigned int n_target_blocks = max_block + 1;
 
92
 
 
93
    std::vector<std::vector<unsigned int> >
 
94
      ndofs(mg_dof.get_tria().n_levels(),
 
95
            std::vector<unsigned int>(n_target_blocks));
 
96
    MGTools::count_dofs_per_block (mg_dof, ndofs, target_component);
 
97
 
 
98
    for (unsigned int level=v.get_minlevel();
 
99
         level<=v.get_maxlevel();++level)
 
100
      {
 
101
        v[level].reinit(n_target_blocks);
 
102
        for (unsigned int b=0; b<n_target_blocks; ++b)
 
103
          v[level].block(b).reinit(ndofs[level][b]);
 
104
        v[level].collect_sizes();
 
105
      }
 
106
  }
 
107
}
 
108
 
 
109
 
 
110
 
 
111
template <class VECTOR>
 
112
template <int dim, class InVector, int spacedim>
 
113
void
 
114
MGTransferPrebuilt<VECTOR>::copy_to_mg (
 
115
  const MGDoFHandler<dim,spacedim>& mg_dof_handler,
 
116
  MGLevelObject<VECTOR>& dst,
 
117
  const InVector& src) const
 
118
{
 
119
  reinit_vector(mg_dof_handler, component_to_block_map, dst);
 
120
  bool first = true;
 
121
  for (unsigned int level=mg_dof_handler.get_tria().n_levels();level != 0;)
 
122
    {
 
123
      --level;
 
124
      VECTOR& dst_level = dst[level];
 
125
 
 
126
      typedef std::vector<std::pair<unsigned int, unsigned int> >::const_iterator IT;
 
127
      for (IT i= copy_indices[level].begin();
 
128
           i != copy_indices[level].end();++i)
 
129
        dst_level(i->second) = src(i->first);
 
130
 
 
131
                                       // For non-DG: degrees of
 
132
                                       // freedom in the refinement
 
133
                                       // face may need special
 
134
                                       // attention, since they belong
 
135
                                       // to the coarse level, but
 
136
                                       // have fine level basis
 
137
                                       // functions
 
138
#ifdef BAERBEL_MG_TEST
 
139
      deallog << "dst [" << level << "] " 
 
140
        << dst[level].l2_norm() << std::endl; 
 
141
#endif
 
142
 
 
143
      if (!first)
 
144
        restrict_and_add (level+1, dst[level], dst[level+1]);
 
145
      first = false;
 
146
#ifdef BAERBEL_MG_TEST
 
147
      deallog << "dst [" << level << "] " 
 
148
        << dst[level].l2_norm() << std::endl; 
 
149
#endif
 
150
    }
 
151
}
 
152
 
 
153
 
 
154
 
 
155
template <typename VECTOR>
34
156
template <int dim, int spacedim>
35
 
void MGTransferPrebuilt<number>::build_matrices (
36
 
  const MGDoFHandler<dim,spacedim> &mg_dof)
 
157
void MGTransferPrebuilt<VECTOR>::build_matrices (
 
158
  const MGDoFHandler<dim,spacedim>           &mg_dof,
 
159
  const std::vector<std::set<unsigned int> > &boundary_indices
 
160
  )
37
161
{
38
162
  const unsigned int n_levels      = mg_dof.get_tria().n_levels();
39
163
  const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
40
 
  
 
164
 
41
165
  sizes.resize(n_levels);
42
166
  for (unsigned int l=0;l<n_levels;++l)
43
167
    sizes[l] = mg_dof.n_dofs(l);
44
 
  
 
168
 
45
169
                                   // reset the size of the array of
46
170
                                   // matrices. call resize(0) first,
47
171
                                   // in order to delete all elements
54
178
                                   // by itself
55
179
  prolongation_matrices.resize (0);
56
180
  prolongation_sparsities.resize (0);
57
 
  
 
181
 
58
182
  for (unsigned int i=0; i<n_levels-1; ++i)
59
183
    {
60
 
      prolongation_sparsities
61
 
        .push_back (std_cxx1x::shared_ptr<SparsityPattern> (new SparsityPattern));
62
 
      prolongation_matrices
63
 
        .push_back (std_cxx1x::shared_ptr<SparseMatrix<double> > (new SparseMatrix<double>));
 
184
      prolongation_sparsities.push_back
 
185
        (std_cxx1x::shared_ptr<SparsityPattern> (new SparsityPattern));
 
186
      prolongation_matrices.push_back
 
187
        (std_cxx1x::shared_ptr<SparseMatrix<double> > (new SparseMatrix<double>));
64
188
    }
65
 
  
 
189
 
66
190
                                   // two fields which will store the
67
191
                                   // indices of the multigrid dofs
68
192
                                   // for a cell and one of its children
69
193
  std::vector<unsigned int> dof_indices_parent (dofs_per_cell);
70
194
  std::vector<unsigned int> dof_indices_child (dofs_per_cell);
71
 
  
 
195
 
72
196
                                   // for each level: first build the sparsity
73
197
                                   // pattern of the matrices and then build the
74
198
                                   // matrices themselves. note that we only
76
200
                                   // level which have children
77
201
  for (unsigned int level=0; level<n_levels-1; ++level)
78
202
    {
 
203
 
79
204
                                       // reset the dimension of the structure.
80
205
                                       // note that for the number of entries
81
206
                                       // per row, the number of parent dofs
89
214
      prolongation_sparsities[level]->reinit (sizes[level+1],
90
215
                                              sizes[level],
91
216
                                              dofs_per_cell+1);
92
 
      
 
217
 
93
218
      for (typename MGDoFHandler<dim,spacedim>::cell_iterator cell=mg_dof.begin(level);
94
219
           cell != mg_dof.end(level); ++cell)
95
220
        if (cell->has_children())
104
229
                                                 // prolongation matrix for
105
230
                                                 // this child
106
231
                const FullMatrix<double> &prolongation
107
 
                  = mg_dof.get_fe().get_prolongation_matrix (child, cell->refinement_case());
 
232
                  = mg_dof.get_fe().get_prolongation_matrix (child,
 
233
                                                             cell->refinement_case());
108
234
 
109
235
                Assert (prolongation.n() != 0, ExcNoProlongation());
110
 
                
 
236
 
111
237
                cell->child(child)->get_mg_dof_indices (dof_indices_child);
112
238
 
113
239
                                                 // now tag the entries in the
116
242
                for (unsigned int i=0; i<dofs_per_cell; ++i)
117
243
                  for (unsigned int j=0; j<dofs_per_cell; ++j)
118
244
                    if (prolongation(i,j) != 0)
119
 
                      {
120
 
                        prolongation_sparsities[level]->add (dof_indices_child[i],
121
 
                                                             dof_indices_parent[j]);
122
 
                      };
123
 
              };
124
 
          };
 
245
                      prolongation_sparsities[level]->add (dof_indices_child[i],
 
246
                                                           dof_indices_parent[j]);
 
247
              }
 
248
          }
 
249
 
125
250
      prolongation_sparsities[level]->compress ();
126
251
 
127
252
      prolongation_matrices[level]->reinit (*prolongation_sparsities[level]);
141
266
                                                 // prolongation matrix for
142
267
                                                 // this child
143
268
                const FullMatrix<double> &prolongation
144
 
                  = mg_dof.get_fe().get_prolongation_matrix (child, cell->refinement_case());
145
 
            
 
269
                  = mg_dof.get_fe().get_prolongation_matrix (child,
 
270
                                                             cell->refinement_case());
 
271
 
146
272
                cell->child(child)->get_mg_dof_indices (dof_indices_child);
147
273
 
148
274
                                                 // now set the entries in the
149
275
                                                 // matrix
150
276
                for (unsigned int i=0; i<dofs_per_cell; ++i)
151
 
                  for (unsigned int j=0; j<dofs_per_cell; ++j)
152
 
                    if (prolongation(i,j) != 0)
153
 
                      prolongation_matrices[level]->set (dof_indices_child[i],
154
 
                                                         dof_indices_parent[j],
155
 
                                                         prolongation(i,j));
 
277
                  prolongation_matrices[level]->set (dof_indices_child[i],
 
278
                                                     dofs_per_cell,
 
279
                                                     &dof_indices_parent[0],
 
280
                                                     &prolongation(i,0),
 
281
                                                     true);
156
282
              }
157
283
          }
158
284
    }
159
285
 
160
 
  copy_indices.resize(mg_dof.get_tria().n_levels());
 
286
 
 
287
                                // impose boundary conditions
 
288
                                // but only in the column of
 
289
                                // the prolongation matrix
 
290
  if (boundary_indices.size() != 0)
 
291
    {
 
292
      std::vector<unsigned int> constrain_indices;
 
293
      for (int level=n_levels-2; level>=0; --level)
 
294
        {
 
295
          if (boundary_indices[level].size() == 0)
 
296
            continue;
 
297
 
 
298
                                // need to delete all the columns in the
 
299
                                // matrix that are on the boundary. to achive
 
300
                                // this, create an array as long as there are
 
301
                                // matrix columns, and find which columns we
 
302
                                // need to filter away.
 
303
          constrain_indices.resize (0);
 
304
          constrain_indices.resize (prolongation_matrices[level]->n(), 0);
 
305
          std::set<unsigned int>::const_iterator dof = boundary_indices[level].begin(),
 
306
            endd = boundary_indices[level].end();
 
307
          for (; dof != endd; ++dof)
 
308
            constrain_indices[*dof] = 1;
 
309
 
 
310
          const unsigned int n_dofs = prolongation_matrices[level]->m();
 
311
          for (unsigned int i=0; i<n_dofs; ++i)
 
312
            {
 
313
              SparseMatrix<double>::iterator
 
314
                start_row = prolongation_matrices[level]->begin(i),
 
315
                end_row   = prolongation_matrices[level]->end(i);
 
316
              for(; start_row != end_row; ++start_row)
 
317
                {
 
318
                  if (constrain_indices[start_row->column()] == 1)
 
319
                    start_row->value() = 0;
 
320
                }
 
321
            }
 
322
        }
 
323
    }
 
324
 
 
325
                                // to find the indices that describe the
 
326
                                // relation between global dofs and local
 
327
                                // numbering on the individual level, first
 
328
                                // create a temp vector where the ith level
 
329
                                // entry contains the respective global
 
330
                                // entry. this gives a neat way to find those
 
331
                                // indices. in a second step, actually build
 
332
                                // the std::vector<std::pair<uint,uint> > that
 
333
                                // only contains the active dofs on the
 
334
                                // levels.
 
335
  interface_dofs.resize(mg_dof.get_tria().n_levels());
 
336
  for(unsigned int l=0; l<mg_dof.get_tria().n_levels(); ++l)
 
337
    interface_dofs[l].resize(mg_dof.n_dofs(l));
 
338
  MGTools::extract_inner_interface_dofs (mg_dof, interface_dofs);
 
339
 
 
340
  copy_indices.resize(n_levels);
 
341
  std::vector<unsigned int> temp_copy_indices;
161
342
  std::vector<unsigned int> global_dof_indices (dofs_per_cell);
162
343
  std::vector<unsigned int> level_dof_indices  (dofs_per_cell);
163
344
  for (int level=mg_dof.get_tria().n_levels()-1; level>=0; --level)
167
348
        level_cell = mg_dof.begin_active(level);
168
349
      const typename MGDoFHandler<dim,spacedim>::active_cell_iterator
169
350
        level_end  = mg_dof.end_active(level);
170
 
      
 
351
 
 
352
      temp_copy_indices.resize (0);
 
353
      temp_copy_indices.resize (mg_dof.n_dofs(level), numbers::invalid_unsigned_int);
 
354
 
171
355
                                       // Compute coarse level right hand side
172
356
                                       // by restricting from fine level.
173
357
      for (; level_cell!=level_end; ++level_cell)
179
363
                                           // numbering
180
364
          global_cell.get_dof_indices(global_dof_indices);
181
365
          level_cell->get_mg_dof_indices (level_dof_indices);
182
 
          
 
366
 
183
367
          for (unsigned int i=0; i<dofs_per_cell; ++i)
184
 
            {
185
 
              copy_indices[level].insert(
186
 
                std::make_pair(global_dof_indices[i], level_dof_indices[i]));
187
 
            }
 
368
          {
 
369
            if(!interface_dofs[level][level_dof_indices[i]])
 
370
              temp_copy_indices[level_dof_indices[i]] = global_dof_indices[i];
 
371
          }
188
372
        }
 
373
 
 
374
                                // now all the active dofs got a valid entry,
 
375
                                // the other ones have an invalid entry. Count
 
376
                                // the invalid entries and then resize the
 
377
                                // copy_indices object. Then, insert the pairs
 
378
                                // of global index and level index into
 
379
                                // copy_indices.
 
380
      const unsigned int n_active_dofs =
 
381
        std::count_if (temp_copy_indices.begin(), temp_copy_indices.end(),
 
382
                       std::bind2nd(std::not_equal_to<unsigned int>(),
 
383
                                    numbers::invalid_unsigned_int));
 
384
      copy_indices[level].resize (n_active_dofs);
 
385
      unsigned int counter = 0;
 
386
      for (unsigned int i=0; i<temp_copy_indices.size(); ++i)
 
387
        if (temp_copy_indices[i] != numbers::invalid_unsigned_int)
 
388
          copy_indices[level][counter++] =
 
389
            std::pair<unsigned int, unsigned int> (temp_copy_indices[i], i);
 
390
      Assert (counter == n_active_dofs, ExcInternalError());
189
391
    }
190
392
}
191
393
 
192
394
 
 
395
//TODO: Use template expander script
193
396
 
194
397
template
195
398
void MGTransferPrebuilt<Vector<float> >::build_matrices<deal_II_dimension>
196
 
(const MGDoFHandler<deal_II_dimension> &mg_dof);
 
399
(const MGDoFHandler<deal_II_dimension> &mg_dof,
 
400
 const std::vector<std::set<unsigned int> >&boundary_indices);
197
401
 
198
402
template
199
403
void MGTransferPrebuilt<Vector<double> >::build_matrices<deal_II_dimension>
200
 
(const MGDoFHandler<deal_II_dimension> &mg_dof);
 
404
(const MGDoFHandler<deal_II_dimension> &mg_dof,
 
405
 const std::vector<std::set<unsigned int> >&boundary_indices);
201
406
 
202
407
template
203
408
void MGTransferPrebuilt<BlockVector<float> >::build_matrices<deal_II_dimension>
204
 
(const MGDoFHandler<deal_II_dimension> &mg_dof);
 
409
(const MGDoFHandler<deal_II_dimension> &mg_dof,
 
410
 const std::vector<std::set<unsigned int> >&boundary_indices);
205
411
 
206
412
template
207
413
void MGTransferPrebuilt<BlockVector<double> >::build_matrices<deal_II_dimension>
208
 
(const MGDoFHandler<deal_II_dimension> &mg_dof);
 
414
(const MGDoFHandler<deal_II_dimension> &mg_dof,
 
415
 const std::vector<std::set<unsigned int> >&boundary_indices);
209
416
 
210
417
template void
211
418
MGTransferPrebuilt<Vector<float> >::copy_to_mg (