24 #ifndef __MESH_FUNCTION_H
25 #define __MESH_FUNCTION_H
31 #include <unordered_set>
32 #include <dolfin/common/Hierarchical.h>
33 #include <dolfin/common/MPI.h>
34 #include <dolfin/common/NoDeleter.h>
35 #include <dolfin/common/Variable.h>
36 #include <dolfin/log/log.h>
37 #include <dolfin/io/File.h>
38 #include "LocalMeshValueCollection.h"
39 #include "MeshDomains.h"
40 #include "MeshEntity.h"
42 #include "MeshConnectivity.h"
97 const std::string filename);
145 std::shared_ptr<const Mesh>
mesh()
const;
293 std::string
str(
bool verbose)
const;
300 std::unique_ptr<T[]> _values;
303 std::shared_ptr<const Mesh> _mesh;
318 template <
typename T>
324 template <
typename T>
326 :
Variable(
"f",
"unnamed MeshFunction"),
332 template <
typename T>
335 :
Variable(
"f",
"unnamed MeshFunction"),
341 template <
typename T>
343 std::size_t dim,
const T& value)
350 template <
typename T>
352 const std::string filename)
353 :
Variable(
"f",
"unnamed MeshFunction"),
356 File file(
mesh->mpi_comm(), filename);
360 template <
typename T>
363 :
Variable(
"f",
"unnamed MeshFunction"),
365 _dim(value_collection.dim()), _size(0)
367 *
this = value_collection;
370 template <
typename T>
373 :
Variable(
"f",
"unnamed MeshFunction"),
376 dolfin_assert(_mesh);
385 set_all(std::numeric_limits<T>::max());
388 const std::size_t D = _mesh->topology().dim();
389 dolfin_assert(
dim <= D);
392 const std::map<std::size_t, std::size_t>& data = domains.
markers(
dim);
395 std::map<std::size_t, std::size_t>::const_iterator it;
396 for (it = data.begin(); it != data.end(); ++it)
399 const std::size_t entity_index = it->first;
400 const T value = it->second;
402 dolfin_assert(entity_index < _size);
403 _values[entity_index] = value;
407 template <
typename T>
409 Variable(
"f",
"unnamed MeshFunction"),
415 template <
typename T>
418 if (_size != f._size)
419 _values.reset(
new T[f._size]);
423 std::copy(f._values.get(), f._values.get() + _size, _values.get());
430 template <
typename T>
433 _dim = mesh_value_collection.
dim();
435 dolfin_assert(_mesh);
438 const std::size_t d = _dim;
439 const std::size_t D = _mesh->topology().dim();
440 dolfin_assert(d <= D);
445 dolfin_assert(!connectivity.
empty());
448 set_all(std::numeric_limits<T>::max());
451 std::unordered_set<std::size_t> entities_values_set;
452 typename std::map<std::pair<std::size_t, std::size_t>, T>::const_iterator it;
453 const std::map<std::pair<std::size_t, std::size_t>, T>& values
454 = mesh_value_collection.
values();
455 for (it = values.begin(); it != values.end(); ++it)
458 const std::size_t cell_index = it->first.first;
459 const std::size_t local_entity = it->first.second;
460 const T value = it->second;
462 std::size_t entity_index = 0;
466 dolfin_assert(cell_index < _mesh->num_cells());
467 entity_index = connectivity(cell_index)[local_entity];
471 entity_index = cell_index;
472 dolfin_assert(local_entity == 0);
476 dolfin_assert(entity_index < _size);
477 _values[entity_index] = value;
480 entities_values_set.insert(entity_index);
484 if (entities_values_set.size() != _size)
485 dolfin_debug(
"Mesh value collection does not contain all values for all entities");
490 template <
typename T>
493 dolfin_assert(_mesh);
497 template <
typename T>
503 template <
typename T>
509 template <
typename T>
515 template <
typename T>
518 return _values.get();
521 template <
typename T>
524 return _values.get();
527 template <
typename T>
530 dolfin_assert(_values);
531 dolfin_assert(&entity.
mesh() == _mesh.get());
532 dolfin_assert(entity.
dim() == _dim);
533 dolfin_assert(entity.
index() < _size);
534 return _values[entity.
index()];
537 template <
typename T>
540 dolfin_assert(_values);
541 dolfin_assert(&entity.
mesh() == _mesh.get());
542 dolfin_assert(entity.
dim() == _dim);
543 dolfin_assert(entity.
index() < _size);
544 return _values[entity.
index()];
547 template <
typename T>
550 dolfin_assert(_values);
551 dolfin_assert(index < _size);
552 return _values[index];
555 template <
typename T>
558 dolfin_assert(_values);
559 dolfin_assert(index < _size);
560 return _values[index];
563 template <
typename T>
571 template <
typename T>
577 "initialize mesh function",
578 "Mesh has not been specified for mesh function");
582 init(_mesh, dim, _mesh->num_entities(dim));
585 template <
typename T>
591 "initialize mesh function",
592 "Mesh has not been specified for mesh function");
595 init(_mesh, dim, size);
598 template <
typename T>
604 init(mesh, dim, mesh->num_entities(dim));
607 template <
typename T>
609 std::size_t dim, std::size_t size)
615 dolfin_assert(mesh->num_entities(dim) == size);
619 _values.reset(
new T[size]);
625 template <
typename T>
628 dolfin_assert(_values);
629 dolfin_assert(index < _size);
630 _values[index] = value;
633 template <
typename T>
636 dolfin_assert(_values);
637 dolfin_assert(_size == values.size());
638 std::copy(values.begin(), values.end(), _values.get());
641 template <
typename T>
646 std::fill(_values.get(), _values.get() + _size, value);
649 template <
typename T>
652 dolfin_assert(_values);
653 std::size_t n = std::count(_values.get(), _values.get() + _size, value);
654 std::vector<std::size_t> indices;
656 for (std::size_t i = 0; i < size(); ++i)
658 if (_values[i] == value)
659 indices.push_back(i);
664 template <
typename T>
670 s << str(
false) << std::endl << std::endl;
671 warning(
"Verbose output of MeshFunctions must be implemented manually.");
682 s <<
"<MeshFunction of topological dimension " << dim()
683 <<
" containing " << size() <<
" values>";