Flow123d
dofhandler.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2007 Technical University of Liberec. All rights reserved.
4  *
5  * Please make a following refer to Flow123d on your project site if you use the program for any purpose,
6  * especially for academic research:
7  * Flow123d, Research Centre: Advanced Remedial Technologies, Technical University of Liberec, Czech Republic
8  *
9  * This program is free software; you can redistribute it and/or modify it under the terms
10  * of the GNU General Public License version 3 as published by the Free Software Foundation.
11  *
12  * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
13  * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
14  * See the GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License along with this program; if not,
17  * write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 021110-1307, USA.
18  *
19  *
20  * $Id$
21  * $Revision$
22  * $LastChangedBy$
23  * $LastChangedDate$
24  *
25  * @file
26  * @brief Declaration of class which handles the ordering of degrees of freedom (dof) and mappings between local and global dofs.
27  * @author Jan Stebel
28  */
29 
30 
31 #include "fem/dofhandler.hh"
32 #include "fem/finite_element.hh"
33 #include "mesh/mesh.h"
34 #include "mesh/partitioning.hh"
35 #include "la/distribution.hh"
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 //template<unsigned int dim, unsigned int spacedim> inline
47 //DOFHandler<dim,spacedim>::DOFHandler(Mesh & _mesh)
48 //: DOFHandlerBase(_mesh),
49 // finite_element(0)
50 //{
51 // object_dofs = new int**[mesh->n_elements()];
52 // for (int i=0; i<mesh->n_elements(); i++)
53 // object_dofs[i] = NULL;
54 //}
55 //
56 //
57 //template<unsigned int dim, unsigned int spacedim> inline
58 //DOFHandler<dim,spacedim>::~DOFHandler()
59 //{
60 // for (int i=0; i<mesh->n_elements(); i++)
61 // if (object_dofs[i] != NULL)
62 // {
63 // for (int j=0; j<mesh->element[i].dim(); j++)
64 // if (object_dofs[i][j] != NULL)
65 // delete[] object_dofs[i][j];
66 //
67 // delete[] object_dofs[i];
68 // }
69 // delete[] object_dofs;
70 //}
71 //
72 //
73 //template<unsigned int dim, unsigned int spacedim> inline
74 //void DOFHandler<dim,spacedim>::distribute_dofs(FiniteElement<dim,spacedim> & fe, const unsigned int offset)
75 //{
76 // // First check if dofs are already distributed.
77 // ASSERT(finite_element == 0, "Attempt to distribute DOFs multiple times!");
78 //
79 // unsigned int next_free_dof = offset;
80 // unsigned int n_obj_dofs[dim+1];
81 //
82 // finite_element = &fe;
83 // global_dof_offset = offset;
84 //
85 // for (unsigned int dm=0; dm <= dim; dm++)
86 // {
87 // n_obj_dofs[dm] = 0;
88 // for (unsigned int m=0; m<dof_multiplicities.size(); m++)
89 // n_obj_dofs[dm] += fe.n_object_dofs(dm, dof_multiplicities[m])*dof_multiplicities[m];
90 // }
91 //
92 // // Broadcast partition of elements to all processes.
93 // int *loc_part;
94 // int myp = mesh->get_part()->get_init_distr()->myp();
95 // if (myp == 0)
96 // {
97 // loc_part = (int*)mesh->get_part()->get_loc_part();
98 // }
99 // else
100 // {
101 // loc_part = new int[mesh->n_elements()];
102 // }
103 // MPI_Bcast(loc_part, mesh->n_elements(), MPI_INT, 0, mesh->get_part()->get_init_distr()->get_comm());
104 //
105 // // Distribute element dofs.
106 // // First we distribute dofs on elements associated to process 0 and so on.
107 // for (int proc=0; proc<mesh->get_part()->get_init_distr()->np(); proc++)
108 // {
109 // if (proc == myp)
110 // loffset_ = next_free_dof;
111 //
112 // FOR_ELEMENTS(mesh, cell)
113 // {
114 // if (cell->dim() != dim) continue;
115 //
116 // if (loc_part[cell.index()] != proc) continue;
117 //
118 // // distribute dofs
119 // // TODO: For the moment we distribute only dofs associated to the cell
120 // // In the future we want to distribute dofs on vertices, lines,
121 // // and triangles as well.
122 // object_dofs[cell.index()] = new int*[dim+1];
123 // for (int i=0; i<dim+1; i++)
124 // object_dofs[cell.index()][i] = NULL;
125 // object_dofs[cell.index()][dim] = new int[n_obj_dofs[dim]];
126 //
127 // for (unsigned int i=0; i<n_obj_dofs[dim]; i++)
128 // object_dofs[cell.index()][dim][i] = next_free_dof++;
129 // }
130 //
131 // if (proc == myp)
132 // lsize_ = next_free_dof - loffset_;
133 // }
134 //
135 // // Finally we free the unused array loc_part.
136 // if (mesh->get_part()->get_init_distr()->myp() != 0)
137 // delete[] loc_part;
138 //
139 //// FOR_ELEMENTS(mesh,cell)
140 //// {
141 //// // skip cells of different dimension
142 //// if (cell->dim() != dim) continue;
143 ////
144 //// // distribute dofs
145 //// // TODO: For the moment we distribute only dofs associated to the cell
146 //// // In the future we want to distribute dofs on vertices, lines,
147 //// // and triangles as well.
148 //// object_dofs[dim][cell] = new int[n_obj_dofs[dim]];
149 //// for (unsigned int i=0; i<n_obj_dofs[dim]; i++)
150 //// object_dofs[dim][cell][i] = next_free_dof++;
151 //// }
152 //
153 //
154 //// FOR_ELEMENTS(mesh,cell)
155 //// {
156 //// // skip cells of different dimension
157 //// if (cell->dim != dim) continue;
158 ////
159 //// // distribute dofs
160 //// for (int dm=0; dm<=dim; dm++)
161 //// {
162 //// for (int i=0; i<cell->n_sides_by_dim(dm); i++)
163 //// {
164 //// void *side = cell->side_by_dim(dm, i);
165 //// // check if node has already assigned dofs, otherwise
166 //// // distribute
167 //// if (object_dofs[dm].find(side) == object_dofs[dm].end())
168 //// {
169 //// object_dofs[dm][side] = new int[n_obj_dofs[dm]];
170 //// for (int i=0; i<n_obj_dofs[dm]; i++)
171 //// object_dofs[dm][side][i] = next_free_dof++;
172 //// }
173 //// }
174 //// }
175 //// }
176 //
177 // n_dofs = next_free_dof - offset;
178 //}
179 //
180 //
181 //
182 //template<unsigned int dim, unsigned int spacedim> inline
183 //const unsigned int DOFHandler<dim,spacedim>::n_local_dofs()
184 //{
185 // return finite_element->n_dofs();
186 //}
187 //
188 //
189 //
190 //template<unsigned int dim, unsigned int spacedim>
191 //void DOFHandler<dim,spacedim>::get_dof_indices(const CellIterator &cell, unsigned int indices[])
192 //{
193 //// void *side;
194 //// unsigned int offset, pid;
195 //
196 // for (unsigned int k=0; k<finite_element->n_object_dofs(dim,DOF_SINGLE); k++)
197 // indices[k] = object_dofs[cell.index()][dim][k];
198 //
199 //// indices.clear();
200 ////
201 //// get_object_dof_indices<0>(cell, indices);
202 //// get_object_dof_indices<1>(cell, indices);
203 //// get_object_dof_indices<2>(cell, indices);
204 //// get_object_dof_indices<3>(cell, indices);
205 //}
206 //
207 ////template<unsigned int dim> template<unsigned int obj_dim> inline void DOFHandler<dim>::get_object_dof_indices(const CellIterator &cell, unsigned int indices[])
208 ////{
209 // // TODO: implement for lower dimensional objects
210 //
211 //// void *side;
212 //// unsigned int offset, pid;
213 ////
214 //// // loop over cell points/lines/triangles/tetrahedra
215 //// for (int i=0; i<n_simplex_objects<dim>(obj_dim); i++)
216 //// {
217 //// side = cell->side_by_dim(obj_dim,i);
218 //// pid = permutation_id<dim,obj_dim>(cell,i);
219 //// offset = 0;
220 //// // loop over dof multiplicities (single dofs, pairs, triples, sextuples)
221 //// for (vector<unsigned int>::iterator m=dof_multiplicities.begin(); m!=dof_multiplicities.end(); m++)
222 //// {
223 //// // loop over particular single dofs/dof pairs/triples/sextuples
224 //// for (int j=0; j<finite_element.n_object_dofs(obj_dim,*m); j++)
225 //// {
226 //// // loop over particular dofs (the single dof/2 dofs in the pair etc.)
227 //// for (int k=0; k<*m; k++)
228 //// indices.push_back(object_dofs[obj_dim][side][offset+Simplex<obj_dim>::pair_permutations[pid][k]]);
229 ////
230 //// offset += *m;
231 //// }
232 //// }
233 //// }
234 ////}
235 //
236 //template<unsigned int dim, unsigned int spacedim> inline
237 //void DOFHandler<dim,spacedim>::get_dof_values(const CellIterator &cell, const Vec &values, double local_values[])
238 //{
239 // unsigned int indices[finite_element->n_dofs()];
240 //
241 // get_dof_indices(cell, indices);
242 // VecGetValues(values, finite_element->n_dofs(), (PetscInt *)indices, local_values);
243 //}
244 
245 
246 
247 DOFHandlerMultiDim::DOFHandlerMultiDim(Mesh& _mesh)
248  : DOFHandlerBase(_mesh),
249  fe1d_(0),
250  fe2d_(0),
251  fe3d_(0)
252 {
253  object_dofs = new int**[mesh_->n_elements()];
254  for (int i=0; i<mesh_->n_elements(); i++)
255  object_dofs[i] = NULL;
256 
257  make_elem_partitioning();
258 }
259 
260 void DOFHandlerMultiDim::distribute_dofs(FiniteElement<1, 3>& fe1d,
261  FiniteElement<2, 3>& fe2d, FiniteElement<3, 3>& fe3d,
262  const unsigned int offset)
263 {
264  // First check if dofs are already distributed.
265  ASSERT((fe1d_ == 0) && (fe2d_ == 0) && (fe3d_ == 0), "Attempt to distribute DOFs multiple times!");
266 
267  unsigned int next_free_dof = offset;
268  // n_obj_dofs[element_dim][general_face_dim]
269  // for every dim of ref. element, number of DoFs on its generalized face (face, edge, vertex)
270  unsigned int n_obj_dofs[4][4];
271 
272 
273  fe1d_ = &fe1d;
274  fe2d_ = &fe2d;
275  fe3d_ = &fe3d;
276  global_dof_offset = offset;
277 
278  for (unsigned int dm=0; dm <= 1; dm++)
279  {
280  n_obj_dofs[1][dm] = 0;
281  for (unsigned int m=0; m<dof_multiplicities.size(); m++)
282  n_obj_dofs[1][dm] += fe1d.n_object_dofs(dm, dof_multiplicities[m])*dof_multiplicities[m];
283  }
284  for (unsigned int dm=0; dm <= 2; dm++)
285  {
286  n_obj_dofs[2][dm] = 0;
287  for (unsigned int m=0; m<dof_multiplicities.size(); m++)
288  n_obj_dofs[2][dm] += fe2d.n_object_dofs(dm, dof_multiplicities[m])*dof_multiplicities[m];
289  }
290  for (unsigned int dm=0; dm <= 3; dm++)
291  {
292  n_obj_dofs[3][dm] = 0;
293  for (unsigned int m=0; m<dof_multiplicities.size(); m++)
294  n_obj_dofs[3][dm] += fe3d.n_object_dofs(dm, dof_multiplicities[m])*dof_multiplicities[m];
295  }
296 
297  // Broadcast partition of elements to all processes.
298  int *loc_part;
299  int myp = mesh_->get_part()->get_init_distr()->myp();
300  if (myp == 0)
301  {
302  loc_part = (int*)mesh_->get_part()->get_loc_part();
303  }
304  else
305  {
306  loc_part = new int[mesh_->n_elements()];
307  }
308  MPI_Bcast(loc_part, mesh_->n_elements(), MPI_INT, 0, mesh_->get_part()->get_init_distr()->get_comm());
309 
310  // Distribute element dofs.
311  // First we distribute dofs on elements associated to process 0 and so on.
312  for (int proc=0; proc<mesh_->get_part()->get_init_distr()->np(); proc++)
313  {
314  if (proc == myp)
315  loffset_ = next_free_dof;
316 
317  FOR_ELEMENTS(mesh_, cell)
318  {
319  if (loc_part[cell.index()] != proc) continue;
320 
321  unsigned int dim = cell->dim();
322 
323  // distribute dofs
324  // TODO: For the moment we distribute only dofs associated to the cell
325  // In the future we want to distribute dofs on vertices, lines,
326  // and triangles as well.
327  object_dofs[cell.index()] = new int*[dim+1];
328  for (int i=0; i<dim+1; i++)
329  object_dofs[cell.index()][i] = NULL;
330  object_dofs[cell.index()][dim] = new int[n_obj_dofs[dim][dim]];
331 
332  for (unsigned int i=0; i<n_obj_dofs[dim][dim]; i++)
333  object_dofs[cell.index()][dim][i] = next_free_dof++;
334  }
335 
336  if (proc == myp) {
337  lsize_ = next_free_dof - loffset_;
338  ds_ = new Distribution(lsize_, PETSC_COMM_WORLD);
339  }
340  }
341 
342  // Finally we free the unused array loc_part.
343  if (mesh_->get_part()->get_init_distr()->myp() != 0)
344  delete[] loc_part;
345 
346  n_dofs = next_free_dof - offset;
347 }
348 
349 void DOFHandlerMultiDim::get_dof_indices(const CellIterator &cell, unsigned int indices[]) const
350 {
351  unsigned int dim = cell->dim();
352  switch (dim)
353  {
354  case 1:
355  for (unsigned int k=0; k<fe1d_->n_object_dofs(dim,DOF_SINGLE); k++)
356  indices[k] = object_dofs[cell.index()][dim][k];
357  break;
358  case 2:
359  for (unsigned int k=0; k<fe2d_->n_object_dofs(dim,DOF_SINGLE); k++)
360  indices[k] = object_dofs[cell.index()][dim][k];
361  break;
362  case 3:
363  for (unsigned int k=0; k<fe3d_->n_object_dofs(dim,DOF_SINGLE); k++)
364  indices[k] = object_dofs[cell.index()][dim][k];
365  break;
366  }
367 }
368 
369 void DOFHandlerMultiDim::get_dof_values(const CellIterator &cell, const Vec &values, double local_values[]) const
370 {
371  int ndofs;
372 
373  switch (cell->dim())
374  {
375  case 1:
376  ndofs = fe1d_->n_dofs();
377  break;
378  case 2:
379  ndofs = fe2d_->n_dofs();
380  break;
381  case 3:
382  ndofs = fe3d_->n_dofs();
383  break;
384  }
385 
386  unsigned int indices[ndofs];
387 
388  get_dof_indices(cell, indices);
389  VecGetValues(values, ndofs, (PetscInt *)indices, local_values);
390 }
391 
392 DOFHandlerMultiDim::~DOFHandlerMultiDim()
393 {
394  for (ElementFullIter elem=mesh_->element.begin(); elem!=mesh_->element.end(); ++elem)
395  if (object_dofs[elem.index()] != NULL)
396  {
397  for (int j=0; j<elem->dim(); j++)
398  if (object_dofs[elem.index()][j] != NULL)
399  delete[] object_dofs[elem.index()][j];
400 
401  delete[] object_dofs[elem.index()];
402  }
403  delete[] object_dofs;
404 }
405 
406 
407 
408 void DOFHandlerMultiDim::make_elem_partitioning()
409 {
410  // create local arrays of elements
411  int *id_4_old = new int[mesh_->n_elements()];
412  int i = 0;
413  FOR_ELEMENTS(mesh_, ele) id_4_old[i++] = ele.index();
414  mesh_->get_part()->id_maps(mesh_->n_elements(), id_4_old, el_ds_, el_4_loc, row_4_el);
415  delete[] id_4_old;
416 
417  // create local array of edges
418  for (unsigned int iedg=0; iedg<mesh_->edges.size(); iedg++)
419  {
420  bool is_edge_local = false;
421  Edge *edg = &mesh_->edges[iedg];
422  for (int sid=0; sid<edg->n_sides; sid++)
423  if (el_is_local(edg->side(sid)->element().index()))
424  {
425  is_edge_local = true;
426  break;
427  }
428  if (is_edge_local)
429  edg_4_loc.push_back(iedg);
430  }
431 
432  // create local array of neighbours
433  for (unsigned int inb=0; inb<mesh_->vb_neighbours_.size(); inb++)
434  {
435  Neighbour *nb = &mesh_->vb_neighbours_[inb];
436  if (el_is_local(mesh_->element.index(nb->element()))
437  || el_is_local(nb->side()->element().index()))
438  nb_4_loc.push_back(inb);
439  }
440 }
441 
442 
443 bool DOFHandlerMultiDim::el_is_local(int index) const
444 {
445  return el_ds_->is_local(row_4_el[index]);
446 }
447 
448 
449 
450 
451 
452 template<> FiniteElement<1,3> *DOFHandlerMultiDim::fe<1>() const { return fe1d_; }
453 template<> FiniteElement<2,3> *DOFHandlerMultiDim::fe<2>() const { return fe2d_; }
454 template<> FiniteElement<3,3> *DOFHandlerMultiDim::fe<3>() const { return fe3d_; }
455 
456 //template class DOFHandler<0,3>;
457 //template class DOFHandler<1,3>;
458 //template class DOFHandler<2,3>;
459 //template class DOFHandler<3,3>;
460 
Generated on Thu May 29 2014 23:14:47 for Flow123d by   doxygen 1.8.4