DOLFIN
DOLFIN C++ interface
MeshFunction.h
1 // Copyright (C) 2006-2009 Anders Logg
2 //
3 // This file is part of DOLFIN.
4 //
5 // DOLFIN is free software: you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published by
7 // the Free Software Foundation, either version 3 of the License, or
8 // (at your option) any later version.
9 //
10 // DOLFIN is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public License
16 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
17 //
18 // Modified by Johan Hoffman, 2007.
19 // Modified by Garth N. Wells, 2010-2013
20 //
21 // First added: 2006-05-22
22 // Last changed: 2013-05-10
23 
24 #ifndef __MESH_FUNCTION_H
25 #define __MESH_FUNCTION_H
26 
27 #include <map>
28 #include <vector>
29 #include <algorithm>
30 
31 #include <memory>
32 #include <unordered_set>
33 #include <dolfin/common/Hierarchical.h>
34 #include <dolfin/common/MPI.h>
35 #include <dolfin/common/NoDeleter.h>
36 #include <dolfin/common/Variable.h>
37 #include <dolfin/log/log.h>
38 #include <dolfin/io/File.h>
39 #include "LocalMeshValueCollection.h"
40 #include "MeshDomains.h"
41 #include "MeshEntity.h"
42 #include "Mesh.h"
43 #include "MeshConnectivity.h"
44 
45 namespace dolfin
46 {
47 
48  class MeshEntity;
49 
56 
57  template <typename T> class MeshFunction : public Variable,
58  public Hierarchical<MeshFunction<T>>
59  {
60  public:
61 
64 
69  explicit MeshFunction(std::shared_ptr<const Mesh> mesh);
70 
77  MeshFunction(std::shared_ptr<const Mesh> mesh, std::size_t dim);
78 
88  MeshFunction(std::shared_ptr<const Mesh> mesh, std::size_t dim,
89  const T& value);
90 
97  MeshFunction(std::shared_ptr<const Mesh> mesh,
98  const std::string filename);
99 
106  MeshFunction(std::shared_ptr<const Mesh> mesh,
107  const MeshValueCollection<T>& value_collection);
108 
117  MeshFunction(std::shared_ptr<const Mesh> mesh,
118  std::size_t dim, const MeshDomains& domains);
119 
125 
128 
135 
141 
146  std::shared_ptr<const Mesh> mesh() const;
147 
152  std::size_t dim() const;
153 
158  bool empty() const;
159 
164  std::size_t size() const;
165 
170  const T* values() const;
171 
176  T* values();
177 
185  T& operator[] (const MeshEntity& entity);
186 
194  const T& operator[] (const MeshEntity& entity) const;
195 
203  T& operator[] (std::size_t index);
204 
212  const T& operator[] (std::size_t index) const;
213 
216  const MeshFunction<T>& operator= (const T& value);
217 
222  void init(std::size_t dim);
223 
231  void init(std::size_t dim, std::size_t size);
232 
239  void init(std::shared_ptr<const Mesh> mesh, std::size_t dim);
240 
250  void init(std::shared_ptr<const Mesh> mesh, std::size_t dim,
251  std::size_t size);
252 
259  void set_value(std::size_t index, const T& value);
260 
262  void set_value(std::size_t index, const T& value, const Mesh& mesh)
263  { set_value(index, value); }
264 
269  void set_values(const std::vector<T>& values);
270 
275  void set_all(const T& value);
276 
285  std::vector<std::size_t> where_equal(T value);
286 
294  std::string str(bool verbose) const;
295 
296  private:
297 
298  // Values at the set of mesh entities. We don't use a
299  // std::vector<T> here because it has trouble with bool, which C++
300  // specialises.
301  std::unique_ptr<T[]> _values;
302 
303  // The mesh
304  std::shared_ptr<const Mesh> _mesh;
305 
306  // Topological dimension
307  std::size_t _dim;
308 
309  // Number of mesh entities
310  std::size_t _size;
311  };
312 
313  template<> std::string MeshFunction<double>::str(bool verbose) const;
314  template<> std::string MeshFunction<std::size_t>::str(bool verbose) const;
315 
316  //---------------------------------------------------------------------------
317  // Implementation of MeshFunction
318  //---------------------------------------------------------------------------
319  template <typename T>
321  {
322  // Do nothing
323  }
324  //---------------------------------------------------------------------------
325  template <typename T>
326  MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh)
327  : Variable("f", "unnamed MeshFunction"),
328  Hierarchical<MeshFunction<T>>(*this), _mesh(mesh), _dim(0), _size(0)
329  {
330  // Do nothing
331  }
332  //---------------------------------------------------------------------------
333  template <typename T>
334  MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh,
335  std::size_t dim)
336  : Variable("f", "unnamed MeshFunction"),
337  Hierarchical<MeshFunction<T>>(*this), _mesh(mesh), _dim(0), _size(0)
338  {
339  init(dim);
340  }
341  //---------------------------------------------------------------------------
342  template <typename T>
343  MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh,
344  std::size_t dim, const T& value)
345  : MeshFunction(mesh, dim)
346 
347  {
348  set_all(value);
349  }
350  //---------------------------------------------------------------------------
351  template <typename T>
352  MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh,
353  const std::string filename)
354  : Variable("f", "unnamed MeshFunction"),
355  Hierarchical<MeshFunction<T>>(*this), _mesh(mesh), _dim(0), _size(0)
356  {
357  File file(mesh->mpi_comm(), filename);
358  file >> *this;
359  }
360  //---------------------------------------------------------------------------
361  template <typename T>
362  MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh,
363  const MeshValueCollection<T>& value_collection)
364  : Variable("f", "unnamed MeshFunction"),
365  Hierarchical<MeshFunction<T>>(*this), _mesh(mesh),
366  _dim(value_collection.dim()), _size(0)
367  {
368  *this = value_collection;
369  }
370  //---------------------------------------------------------------------------
371  template <typename T>
372  MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh,
373  std::size_t dim, const MeshDomains& domains)
374  : Variable("f", "unnamed MeshFunction"),
375  Hierarchical<MeshFunction<T>>(*this), _mesh(mesh), _dim(0), _size(0)
376  {
377  dolfin_assert(_mesh);
378 
379  // Initialise MeshFunction
380  init(dim);
381 
382  // Initialise mesh
383  mesh->init(dim);
384 
385  // Set MeshFunction with default value
386  set_all(std::numeric_limits<T>::max());
387 
388  // Get mesh dimension
389  const std::size_t D = _mesh->topology().dim();
390  dolfin_assert(dim <= D);
391 
392  // Get domain data
393  const std::map<std::size_t, std::size_t>& data = domains.markers(dim);
394 
395  // Iterate over all values and copy into MeshFunctions
396  std::map<std::size_t, std::size_t>::const_iterator it;
397  for (it = data.begin(); it != data.end(); ++it)
398  {
399  // Get value collection entry data
400  const std::size_t entity_index = it->first;
401  const T value = it->second;
402 
403  dolfin_assert(entity_index < _size);
404  _values[entity_index] = value;
405  }
406  }
407  //---------------------------------------------------------------------------
408  template <typename T>
410  Variable("f", "unnamed MeshFunction"),
411  Hierarchical<MeshFunction<T>>(*this), _dim(0), _size(0)
412  {
413  *this = f;
414  }
415  //---------------------------------------------------------------------------
416  template <typename T>
418  {
419  if (_size != f._size)
420  _values.reset(new T[f._size]);
421  _mesh = f._mesh;
422  _dim = f._dim;
423  _size = f._size;
424  std::copy(f._values.get(), f._values.get() + _size, _values.get());
425 
426  Hierarchical<MeshFunction<T>>::operator=(f);
427 
428  return *this;
429  }
430  //---------------------------------------------------------------------------
431  template <typename T>
433  {
434  _dim = mesh_value_collection.dim();
435  init(_dim);
436  dolfin_assert(_mesh);
437 
438  // Get mesh connectivity D --> d
439  const std::size_t d = _dim;
440  const std::size_t D = _mesh->topology().dim();
441  dolfin_assert(d <= D);
442 
443  // Generate connectivity if it does not exist
444  _mesh->init(D, d);
445  const MeshConnectivity& connectivity = _mesh->topology()(D, d);
446  dolfin_assert(!connectivity.empty());
447 
448  // Set MeshFunction with default value
449  set_all(std::numeric_limits<T>::max());
450 
451  // Iterate over all values
452  std::unordered_set<std::size_t> entities_values_set;
453  typename std::map<std::pair<std::size_t, std::size_t>, T>::const_iterator it;
454  const std::map<std::pair<std::size_t, std::size_t>, T>& values
455  = mesh_value_collection.values();
456  for (it = values.begin(); it != values.end(); ++it)
457  {
458  // Get value collection entry data
459  const std::size_t cell_index = it->first.first;
460  const std::size_t local_entity = it->first.second;
461  const T value = it->second;
462 
463  std::size_t entity_index = 0;
464  if (d != D)
465  {
466  // Get global (local to to process) entity index
467  dolfin_assert(cell_index < _mesh->num_cells());
468  entity_index = connectivity(cell_index)[local_entity];
469  }
470  else
471  {
472  entity_index = cell_index;
473  dolfin_assert(local_entity == 0);
474  }
475 
476  // Set value for entity
477  dolfin_assert(entity_index < _size);
478  _values[entity_index] = value;
479 
480  // Add entity index to set (used to check that all values are set)
481  entities_values_set.insert(entity_index);
482  }
483 
484  // Check that all values have been set, if not issue a debug message
485  if (entities_values_set.size() != _size)
486  dolfin_debug("Mesh value collection does not contain all values for all entities");
487 
488  return *this;
489  }
490  //---------------------------------------------------------------------------
491  template <typename T>
492  std::shared_ptr<const Mesh> MeshFunction<T>::mesh() const
493  {
494  dolfin_assert(_mesh);
495  return _mesh;
496  }
497  //---------------------------------------------------------------------------
498  template <typename T>
499  std::size_t MeshFunction<T>::dim() const
500  {
501  return _dim;
502  }
503  //---------------------------------------------------------------------------
504  template <typename T>
506  {
507  return _size == 0;
508  }
509  //---------------------------------------------------------------------------
510  template <typename T>
511  std::size_t MeshFunction<T>::size() const
512  {
513  return _size;
514  }
515  //---------------------------------------------------------------------------
516  template <typename T>
517  const T* MeshFunction<T>::values() const
518  {
519  return _values.get();
520  }
521  //---------------------------------------------------------------------------
522  template <typename T>
524  {
525  return _values.get();
526  }
527  //---------------------------------------------------------------------------
528  template <typename T>
530  {
531  dolfin_assert(_values);
532  dolfin_assert(&entity.mesh() == _mesh.get());
533  dolfin_assert(entity.dim() == _dim);
534  dolfin_assert(entity.index() < _size);
535  return _values[entity.index()];
536  }
537  //---------------------------------------------------------------------------
538  template <typename T>
539  const T& MeshFunction<T>::operator[] (const MeshEntity& entity) const
540  {
541  dolfin_assert(_values);
542  dolfin_assert(&entity.mesh() == _mesh.get());
543  dolfin_assert(entity.dim() == _dim);
544  dolfin_assert(entity.index() < _size);
545  return _values[entity.index()];
546  }
547  //---------------------------------------------------------------------------
548  template <typename T>
549  T& MeshFunction<T>::operator[] (std::size_t index)
550  {
551  dolfin_assert(_values);
552  dolfin_assert(index < _size);
553  return _values[index];
554  }
555  //---------------------------------------------------------------------------
556  template <typename T>
557  const T& MeshFunction<T>::operator[] (std::size_t index) const
558  {
559  dolfin_assert(_values);
560  dolfin_assert(index < _size);
561  return _values[index];
562  }
563  //---------------------------------------------------------------------------
564  template <typename T>
566  {
567  set_all(value);
568  //Hierarchical<MeshFunction<T>>::operator=(value);
569  return *this;
570  }
571  //---------------------------------------------------------------------------
572  template <typename T>
573  void MeshFunction<T>::init(std::size_t dim)
574  {
575  if (!_mesh)
576  {
577  dolfin_error("MeshFunction.h",
578  "initialize mesh function",
579  "Mesh has not been specified for mesh function");
580 
581  }
582  _mesh->init(dim);
583  init(_mesh, dim, _mesh->num_entities(dim));
584  }
585  //---------------------------------------------------------------------------
586  template <typename T>
587  void MeshFunction<T>::init(std::size_t dim, std::size_t size)
588  {
589  if (!_mesh)
590  {
591  dolfin_error("MeshFunction.h",
592  "initialize mesh function",
593  "Mesh has not been specified for mesh function");
594  }
595  _mesh->init(dim);
596  init(_mesh, dim, size);
597  }
598  //---------------------------------------------------------------------------
599  template <typename T>
600  void MeshFunction<T>::init(std::shared_ptr<const Mesh> mesh,
601  std::size_t dim)
602  {
603  dolfin_assert(mesh);
604  mesh->init(dim);
605  init(mesh, dim, mesh->num_entities(dim));
606  }
607  //---------------------------------------------------------------------------
608  template <typename T>
609  void MeshFunction<T>::init(std::shared_ptr<const Mesh> mesh,
610  std::size_t dim, std::size_t size)
611  {
612  dolfin_assert(mesh);
613 
614  // Initialize mesh for entities of given dimension
615  mesh->init(dim);
616  dolfin_assert(mesh->num_entities(dim) == size);
617 
618  // Initialize data
619  if (_size != size)
620  _values.reset(new T[size]);
621  _mesh = mesh;
622  _dim = dim;
623  _size = size;
624  }
625  //---------------------------------------------------------------------------
626  template <typename T>
627  void MeshFunction<T>::set_value(std::size_t index, const T& value)
628  {
629  dolfin_assert(_values);
630  dolfin_assert(index < _size);
631  _values[index] = value;
632  }
633  //---------------------------------------------------------------------------
634  template <typename T>
635  void MeshFunction<T>::set_values(const std::vector<T>& values)
636  {
637  dolfin_assert(_values);
638  dolfin_assert(_size == values.size());
639  std::copy(values.begin(), values.end(), _values.get());
640  }
641  //---------------------------------------------------------------------------
642  template <typename T>
643  void MeshFunction<T>::set_all(const T& value)
644  {
645  //dolfin_assert(_values);
646  if(_values)
647  std::fill(_values.get(), _values.get() + _size, value);
648  }
649  //---------------------------------------------------------------------------
650  template <typename T>
651  std::vector<std::size_t> MeshFunction<T>::where_equal(T value)
652  {
653  dolfin_assert(_values);
654  std::size_t n = std::count(_values.get(), _values.get() + _size, value);
655  std::vector<std::size_t> indices;
656  indices.reserve(n);
657  for (std::size_t i = 0; i < size(); ++i)
658  {
659  if (_values[i] == value)
660  indices.push_back(i);
661  }
662  return indices;
663  }
664  //---------------------------------------------------------------------------
665  template <typename T>
666  std::string MeshFunction<T>::str(bool verbose) const
667  {
668  std::stringstream s;
669  if (verbose)
670  {
671  s << str(false) << std::endl << std::endl;
672  warning("Verbose output of MeshFunctions must be implemented manually.");
673 
674  // This has been disabled as it severely restricts the ease with which
675  // templated MeshFunctions can be used, e.g. it is not possible to
676  // template over std::vector.
677 
678  //for (std::size_t i = 0; i < _size; i++)
679  // s << " (" << _dim << ", " << i << "): " << _values[i] << std::endl;
680  }
681  else
682  {
683  s << "<MeshFunction of topological dimension " << dim()
684  << " containing " << size() << " values>";
685  }
686 
687  return s.str();
688  }
689  //---------------------------------------------------------------------------
690 
691 }
692 
693 #endif
Definition: File.h:46
Definition: Hierarchical.h:44
Definition: MeshConnectivity.h:40
bool empty() const
Return true if the total number of connections is equal to zero.
Definition: MeshConnectivity.h:56
Definition: MeshDomains.h:42
std::map< std::size_t, std::size_t > & markers(std::size_t dim)
Definition: MeshDomains.cpp:62
Definition: MeshEntity.h:43
std::size_t dim() const
Definition: MeshEntity.h:106
std::size_t index() const
Definition: MeshEntity.h:113
const Mesh & mesh() const
Definition: MeshEntity.h:99
Definition: MeshFunction.h:59
MeshFunction< T > & operator=(const MeshValueCollection< T > &mesh)
Definition: MeshFunction.h:432
std::shared_ptr< const Mesh > mesh() const
Definition: MeshFunction.h:492
std::size_t dim() const
Definition: MeshFunction.h:499
MeshFunction(std::shared_ptr< const Mesh > mesh, std::size_t dim, const MeshDomains &domains)
Definition: MeshFunction.h:372
void set_all(const T &value)
Definition: MeshFunction.h:643
MeshFunction(const MeshFunction< T > &f)
Definition: MeshFunction.h:409
bool empty() const
Definition: MeshFunction.h:505
~MeshFunction()
Destructor.
Definition: MeshFunction.h:127
T * values()
Definition: MeshFunction.h:523
void init(std::shared_ptr< const Mesh > mesh, std::size_t dim, std::size_t size)
Definition: MeshFunction.h:609
MeshFunction(std::shared_ptr< const Mesh > mesh, const std::string filename)
Definition: MeshFunction.h:352
std::size_t size() const
Definition: MeshFunction.h:511
std::string str(bool verbose) const
Definition: MeshFunction.h:666
MeshFunction(std::shared_ptr< const Mesh > mesh, std::size_t dim, const T &value)
Definition: MeshFunction.h:343
void init(std::size_t dim)
Definition: MeshFunction.h:573
MeshFunction()
Create empty mesh function.
Definition: MeshFunction.h:320
MeshFunction< T > & operator=(const MeshFunction< T > &f)
Definition: MeshFunction.h:417
void set_value(std::size_t index, const T &value)
Definition: MeshFunction.h:627
T & operator[](const MeshEntity &entity)
Definition: MeshFunction.h:529
MeshFunction(std::shared_ptr< const Mesh > mesh)
Definition: MeshFunction.h:326
void set_values(const std::vector< T > &values)
Definition: MeshFunction.h:635
void init(std::size_t dim, std::size_t size)
Definition: MeshFunction.h:587
MeshFunction(std::shared_ptr< const Mesh > mesh, const MeshValueCollection< T > &value_collection)
Definition: MeshFunction.h:362
void set_value(std::size_t index, const T &value, const Mesh &mesh)
Compatibility function for use in SubDomains.
Definition: MeshFunction.h:262
std::vector< std::size_t > where_equal(T value)
Definition: MeshFunction.h:651
void init(std::shared_ptr< const Mesh > mesh, std::size_t dim)
Definition: MeshFunction.h:600
const T * values() const
Definition: MeshFunction.h:517
MeshFunction(std::shared_ptr< const Mesh > mesh, std::size_t dim)
Definition: MeshFunction.h:334
Definition: MeshValueCollection.h:51
std::size_t dim() const
Definition: MeshValueCollection.h:389
std::map< std::pair< std::size_t, std::size_t >, T > & values()
Definition: MeshValueCollection.h:521
Definition: Mesh.h:84
Common base class for DOLFIN variables.
Definition: Variable.h:36
Definition: adapt.h:30
void warning(std::string msg,...)
Print warning.
Definition: log.cpp:115
void dolfin_error(std::string location, std::string task, std::string reason,...)
Definition: log.cpp:129
void init(int argc, char *argv[])
Definition: init.cpp:27