24 #ifndef __MESH_FUNCTION_H
25 #define __MESH_FUNCTION_H
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"
43 #include "MeshConnectivity.h"
98 const std::string filename);
146 std::shared_ptr<const Mesh>
mesh()
const;
294 std::string
str(
bool verbose)
const;
301 std::unique_ptr<T[]> _values;
304 std::shared_ptr<const Mesh> _mesh;
319 template <
typename T>
325 template <
typename T>
327 :
Variable(
"f",
"unnamed MeshFunction"),
333 template <
typename T>
336 :
Variable(
"f",
"unnamed MeshFunction"),
342 template <
typename T>
344 std::size_t dim,
const T& value)
351 template <
typename T>
353 const std::string filename)
354 :
Variable(
"f",
"unnamed MeshFunction"),
357 File file(
mesh->mpi_comm(), filename);
361 template <
typename T>
364 :
Variable(
"f",
"unnamed MeshFunction"),
366 _dim(value_collection.dim()), _size(0)
368 *
this = value_collection;
371 template <
typename T>
374 :
Variable(
"f",
"unnamed MeshFunction"),
377 dolfin_assert(_mesh);
386 set_all(std::numeric_limits<T>::max());
389 const std::size_t D = _mesh->topology().dim();
390 dolfin_assert(
dim <= D);
393 const std::map<std::size_t, std::size_t>& data = domains.
markers(
dim);
396 std::map<std::size_t, std::size_t>::const_iterator it;
397 for (it = data.begin(); it != data.end(); ++it)
400 const std::size_t entity_index = it->first;
401 const T value = it->second;
403 dolfin_assert(entity_index < _size);
404 _values[entity_index] = value;
408 template <
typename T>
410 Variable(
"f",
"unnamed MeshFunction"),
416 template <
typename T>
419 if (_size != f._size)
420 _values.reset(
new T[f._size]);
424 std::copy(f._values.get(), f._values.get() + _size, _values.get());
431 template <
typename T>
434 _dim = mesh_value_collection.
dim();
436 dolfin_assert(_mesh);
439 const std::size_t d = _dim;
440 const std::size_t D = _mesh->topology().dim();
441 dolfin_assert(d <= D);
446 dolfin_assert(!connectivity.
empty());
449 set_all(std::numeric_limits<T>::max());
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)
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;
463 std::size_t entity_index = 0;
467 dolfin_assert(cell_index < _mesh->num_cells());
468 entity_index = connectivity(cell_index)[local_entity];
472 entity_index = cell_index;
473 dolfin_assert(local_entity == 0);
477 dolfin_assert(entity_index < _size);
478 _values[entity_index] = value;
481 entities_values_set.insert(entity_index);
485 if (entities_values_set.size() != _size)
486 dolfin_debug(
"Mesh value collection does not contain all values for all entities");
491 template <
typename T>
494 dolfin_assert(_mesh);
498 template <
typename T>
504 template <
typename T>
510 template <
typename T>
516 template <
typename T>
519 return _values.get();
522 template <
typename T>
525 return _values.get();
528 template <
typename T>
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()];
538 template <
typename T>
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()];
548 template <
typename T>
551 dolfin_assert(_values);
552 dolfin_assert(index < _size);
553 return _values[index];
556 template <
typename T>
559 dolfin_assert(_values);
560 dolfin_assert(index < _size);
561 return _values[index];
564 template <
typename T>
572 template <
typename T>
578 "initialize mesh function",
579 "Mesh has not been specified for mesh function");
583 init(_mesh, dim, _mesh->num_entities(dim));
586 template <
typename T>
592 "initialize mesh function",
593 "Mesh has not been specified for mesh function");
596 init(_mesh, dim, size);
599 template <
typename T>
605 init(mesh, dim, mesh->num_entities(dim));
608 template <
typename T>
610 std::size_t dim, std::size_t size)
616 dolfin_assert(mesh->num_entities(dim) == size);
620 _values.reset(
new T[size]);
626 template <
typename T>
629 dolfin_assert(_values);
630 dolfin_assert(index < _size);
631 _values[index] = value;
634 template <
typename T>
637 dolfin_assert(_values);
638 dolfin_assert(_size == values.size());
639 std::copy(values.begin(), values.end(), _values.get());
642 template <
typename T>
647 std::fill(_values.get(), _values.get() + _size, value);
650 template <
typename T>
653 dolfin_assert(_values);
654 std::size_t n = std::count(_values.get(), _values.get() + _size, value);
655 std::vector<std::size_t> indices;
657 for (std::size_t i = 0; i < size(); ++i)
659 if (_values[i] == value)
660 indices.push_back(i);
665 template <
typename T>
671 s << str(
false) << std::endl << std::endl;
672 warning(
"Verbose output of MeshFunctions must be implemented manually.");
683 s <<
"<MeshFunction of topological dimension " << dim()
684 <<
" containing " << size() <<
" values>";