DOLFIN
DOLFIN C++ interface
Loading...
Searching...
No Matches
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#include <memory>
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"
41#include "Mesh.h"
42#include "MeshConnectivity.h"
43
44namespace dolfin
45{
46
47 class MeshEntity;
48
55
56 template <typename T> class MeshFunction : public Variable,
57 public Hierarchical<MeshFunction<T>>
58 {
59 public:
60
63
68 explicit MeshFunction(std::shared_ptr<const Mesh> mesh);
69
76 MeshFunction(std::shared_ptr<const Mesh> mesh, std::size_t dim);
77
87 MeshFunction(std::shared_ptr<const Mesh> mesh, std::size_t dim,
88 const T& value);
89
96 MeshFunction(std::shared_ptr<const Mesh> mesh,
97 const std::string filename);
98
105 MeshFunction(std::shared_ptr<const Mesh> mesh,
106 const MeshValueCollection<T>& value_collection);
107
116 MeshFunction(std::shared_ptr<const Mesh> mesh,
117 std::size_t dim, const MeshDomains& domains);
118
124
127
134
140
145 std::shared_ptr<const Mesh> mesh() const;
146
151 std::size_t dim() const;
152
157 bool empty() const;
158
163 std::size_t size() const;
164
169 const T* values() const;
170
175 T* values();
176
184 T& operator[] (const MeshEntity& entity);
185
193 const T& operator[] (const MeshEntity& entity) const;
194
202 T& operator[] (std::size_t index);
203
211 const T& operator[] (std::size_t index) const;
212
215 const MeshFunction<T>& operator= (const T& value);
216
221 void init(std::size_t dim);
222
230 void init(std::size_t dim, std::size_t size);
231
238 void init(std::shared_ptr<const Mesh> mesh, std::size_t dim);
239
249 void init(std::shared_ptr<const Mesh> mesh, std::size_t dim,
250 std::size_t size);
251
258 void set_value(std::size_t index, const T& value);
259
261 void set_value(std::size_t index, const T& value, const Mesh& mesh)
262 { set_value(index, value); }
263
268 void set_values(const std::vector<T>& values);
269
274 void set_all(const T& value);
275
284 std::vector<std::size_t> where_equal(T value);
285
293 std::string str(bool verbose) const;
294
295 private:
296
297 // Values at the set of mesh entities. We don't use a
298 // std::vector<T> here because it has trouble with bool, which C++
299 // specialises.
300 std::unique_ptr<T[]> _values;
301
302 // The mesh
303 std::shared_ptr<const Mesh> _mesh;
304
305 // Topological dimension
306 std::size_t _dim;
307
308 // Number of mesh entities
309 std::size_t _size;
310 };
311
312 template<> std::string MeshFunction<double>::str(bool verbose) const;
313 template<> std::string MeshFunction<std::size_t>::str(bool verbose) const;
314
315 //---------------------------------------------------------------------------
316 // Implementation of MeshFunction
317 //---------------------------------------------------------------------------
318 template <typename T>
320 {
321 // Do nothing
322 }
323 //---------------------------------------------------------------------------
324 template <typename T>
325 MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh)
326 : Variable("f", "unnamed MeshFunction"),
327 Hierarchical<MeshFunction<T>>(*this), _mesh(mesh), _dim(0), _size(0)
328 {
329 // Do nothing
330 }
331 //---------------------------------------------------------------------------
332 template <typename T>
333 MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh,
334 std::size_t dim)
335 : Variable("f", "unnamed MeshFunction"),
336 Hierarchical<MeshFunction<T>>(*this), _mesh(mesh), _dim(0), _size(0)
337 {
338 init(dim);
339 }
340 //---------------------------------------------------------------------------
341 template <typename T>
342 MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh,
343 std::size_t dim, const T& value)
344 : MeshFunction(mesh, dim)
345
346 {
347 set_all(value);
348 }
349 //---------------------------------------------------------------------------
350 template <typename T>
351 MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh,
352 const std::string filename)
353 : Variable("f", "unnamed MeshFunction"),
354 Hierarchical<MeshFunction<T>>(*this), _mesh(mesh), _dim(0), _size(0)
355 {
356 File file(mesh->mpi_comm(), filename);
357 file >> *this;
358 }
359 //---------------------------------------------------------------------------
360 template <typename T>
361 MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh,
362 const MeshValueCollection<T>& value_collection)
363 : Variable("f", "unnamed MeshFunction"),
364 Hierarchical<MeshFunction<T>>(*this), _mesh(mesh),
365 _dim(value_collection.dim()), _size(0)
366 {
367 *this = value_collection;
368 }
369 //---------------------------------------------------------------------------
370 template <typename T>
371 MeshFunction<T>::MeshFunction(std::shared_ptr<const Mesh> mesh,
372 std::size_t dim, const MeshDomains& domains)
373 : Variable("f", "unnamed MeshFunction"),
374 Hierarchical<MeshFunction<T>>(*this), _mesh(mesh), _dim(0), _size(0)
375 {
376 dolfin_assert(_mesh);
377
378 // Initialise MeshFunction
379 init(dim);
380
381 // Initialise mesh
382 mesh->init(dim);
383
384 // Set MeshFunction with default value
385 set_all(std::numeric_limits<T>::max());
386
387 // Get mesh dimension
388 const std::size_t D = _mesh->topology().dim();
389 dolfin_assert(dim <= D);
390
391 // Get domain data
392 const std::map<std::size_t, std::size_t>& data = domains.markers(dim);
393
394 // Iterate over all values and copy into MeshFunctions
395 std::map<std::size_t, std::size_t>::const_iterator it;
396 for (it = data.begin(); it != data.end(); ++it)
397 {
398 // Get value collection entry data
399 const std::size_t entity_index = it->first;
400 const T value = it->second;
401
402 dolfin_assert(entity_index < _size);
403 _values[entity_index] = value;
404 }
405 }
406 //---------------------------------------------------------------------------
407 template <typename T>
409 Variable("f", "unnamed MeshFunction"),
410 Hierarchical<MeshFunction<T>>(*this), _dim(0), _size(0)
411 {
412 *this = f;
413 }
414 //---------------------------------------------------------------------------
415 template <typename T>
417 {
418 if (_size != f._size)
419 _values.reset(new T[f._size]);
420 _mesh = f._mesh;
421 _dim = f._dim;
422 _size = f._size;
423 std::copy(f._values.get(), f._values.get() + _size, _values.get());
424
425 Hierarchical<MeshFunction<T>>::operator=(f);
426
427 return *this;
428 }
429 //---------------------------------------------------------------------------
430 template <typename T>
432 {
433 _dim = mesh_value_collection.dim();
434 init(_dim);
435 dolfin_assert(_mesh);
436
437 // Get mesh connectivity D --> d
438 const std::size_t d = _dim;
439 const std::size_t D = _mesh->topology().dim();
440 dolfin_assert(d <= D);
441
442 // Generate connectivity if it does not exist
443 _mesh->init(D, d);
444 const MeshConnectivity& connectivity = _mesh->topology()(D, d);
445 dolfin_assert(!connectivity.empty());
446
447 // Set MeshFunction with default value
448 set_all(std::numeric_limits<T>::max());
449
450 // Iterate over all values
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)
456 {
457 // Get value collection entry data
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;
461
462 std::size_t entity_index = 0;
463 if (d != D)
464 {
465 // Get global (local to to process) entity index
466 dolfin_assert(cell_index < _mesh->num_cells());
467 entity_index = connectivity(cell_index)[local_entity];
468 }
469 else
470 {
471 entity_index = cell_index;
472 dolfin_assert(local_entity == 0);
473 }
474
475 // Set value for entity
476 dolfin_assert(entity_index < _size);
477 _values[entity_index] = value;
478
479 // Add entity index to set (used to check that all values are set)
480 entities_values_set.insert(entity_index);
481 }
482
483 // Check that all values have been set, if not issue a debug message
484 if (entities_values_set.size() != _size)
485 dolfin_debug("Mesh value collection does not contain all values for all entities");
486
487 return *this;
488 }
489 //---------------------------------------------------------------------------
490 template <typename T>
491 std::shared_ptr<const Mesh> MeshFunction<T>::mesh() const
492 {
493 dolfin_assert(_mesh);
494 return _mesh;
495 }
496 //---------------------------------------------------------------------------
497 template <typename T>
498 std::size_t MeshFunction<T>::dim() const
499 {
500 return _dim;
501 }
502 //---------------------------------------------------------------------------
503 template <typename T>
505 {
506 return _size == 0;
507 }
508 //---------------------------------------------------------------------------
509 template <typename T>
510 std::size_t MeshFunction<T>::size() const
511 {
512 return _size;
513 }
514 //---------------------------------------------------------------------------
515 template <typename T>
516 const T* MeshFunction<T>::values() const
517 {
518 return _values.get();
519 }
520 //---------------------------------------------------------------------------
521 template <typename T>
523 {
524 return _values.get();
525 }
526 //---------------------------------------------------------------------------
527 template <typename T>
529 {
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()];
535 }
536 //---------------------------------------------------------------------------
537 template <typename T>
538 const T& MeshFunction<T>::operator[] (const MeshEntity& entity) const
539 {
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()];
545 }
546 //---------------------------------------------------------------------------
547 template <typename T>
548 T& MeshFunction<T>::operator[] (std::size_t index)
549 {
550 dolfin_assert(_values);
551 dolfin_assert(index < _size);
552 return _values[index];
553 }
554 //---------------------------------------------------------------------------
555 template <typename T>
556 const T& MeshFunction<T>::operator[] (std::size_t index) const
557 {
558 dolfin_assert(_values);
559 dolfin_assert(index < _size);
560 return _values[index];
561 }
562 //---------------------------------------------------------------------------
563 template <typename T>
565 {
566 set_all(value);
567 //Hierarchical<MeshFunction<T>>::operator=(value);
568 return *this;
569 }
570 //---------------------------------------------------------------------------
571 template <typename T>
572 void MeshFunction<T>::init(std::size_t dim)
573 {
574 if (!_mesh)
575 {
576 dolfin_error("MeshFunction.h",
577 "initialize mesh function",
578 "Mesh has not been specified for mesh function");
579
580 }
581 _mesh->init(dim);
582 init(_mesh, dim, _mesh->num_entities(dim));
583 }
584 //---------------------------------------------------------------------------
585 template <typename T>
586 void MeshFunction<T>::init(std::size_t dim, std::size_t size)
587 {
588 if (!_mesh)
589 {
590 dolfin_error("MeshFunction.h",
591 "initialize mesh function",
592 "Mesh has not been specified for mesh function");
593 }
594 _mesh->init(dim);
595 init(_mesh, dim, size);
596 }
597 //---------------------------------------------------------------------------
598 template <typename T>
599 void MeshFunction<T>::init(std::shared_ptr<const Mesh> mesh,
600 std::size_t dim)
601 {
602 dolfin_assert(mesh);
603 mesh->init(dim);
604 init(mesh, dim, mesh->num_entities(dim));
605 }
606 //---------------------------------------------------------------------------
607 template <typename T>
608 void MeshFunction<T>::init(std::shared_ptr<const Mesh> mesh,
609 std::size_t dim, std::size_t size)
610 {
611 dolfin_assert(mesh);
612
613 // Initialize mesh for entities of given dimension
614 mesh->init(dim);
615 dolfin_assert(mesh->num_entities(dim) == size);
616
617 // Initialize data
618 if (_size != size)
619 _values.reset(new T[size]);
620 _mesh = mesh;
621 _dim = dim;
622 _size = size;
623 }
624 //---------------------------------------------------------------------------
625 template <typename T>
626 void MeshFunction<T>::set_value(std::size_t index, const T& value)
627 {
628 dolfin_assert(_values);
629 dolfin_assert(index < _size);
630 _values[index] = value;
631 }
632 //---------------------------------------------------------------------------
633 template <typename T>
634 void MeshFunction<T>::set_values(const std::vector<T>& values)
635 {
636 dolfin_assert(_values);
637 dolfin_assert(_size == values.size());
638 std::copy(values.begin(), values.end(), _values.get());
639 }
640 //---------------------------------------------------------------------------
641 template <typename T>
642 void MeshFunction<T>::set_all(const T& value)
643 {
644 //dolfin_assert(_values);
645 if(_values)
646 std::fill(_values.get(), _values.get() + _size, value);
647 }
648 //---------------------------------------------------------------------------
649 template <typename T>
650 std::vector<std::size_t> MeshFunction<T>::where_equal(T value)
651 {
652 dolfin_assert(_values);
653 std::size_t n = std::count(_values.get(), _values.get() + _size, value);
654 std::vector<std::size_t> indices;
655 indices.reserve(n);
656 for (std::size_t i = 0; i < size(); ++i)
657 {
658 if (_values[i] == value)
659 indices.push_back(i);
660 }
661 return indices;
662 }
663 //---------------------------------------------------------------------------
664 template <typename T>
665 std::string MeshFunction<T>::str(bool verbose) const
666 {
667 std::stringstream s;
668 if (verbose)
669 {
670 s << str(false) << std::endl << std::endl;
671 warning("Verbose output of MeshFunctions must be implemented manually.");
672
673 // This has been disabled as it severely restricts the ease with which
674 // templated MeshFunctions can be used, e.g. it is not possible to
675 // template over std::vector.
676
677 //for (std::size_t i = 0; i < _size; i++)
678 // s << " (" << _dim << ", " << i << "): " << _values[i] << std::endl;
679 }
680 else
681 {
682 s << "<MeshFunction of topological dimension " << dim()
683 << " containing " << size() << " values>";
684 }
685
686 return s.str();
687 }
688 //---------------------------------------------------------------------------
689
690}
691
692#endif
Definition File.h:46
Definition Hierarchical.h:44
Definition MeshConnectivity.h:42
bool empty() const
Return true if the total number of connections is equal to zero.
Definition MeshConnectivity.h:58
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
const Mesh & mesh() const
Definition MeshEntity.h:99
std::size_t dim() const
Definition MeshEntity.h:106
std::size_t index() const
Definition MeshEntity.h:113
Definition MeshFunction.h:58
MeshFunction< T > & operator=(const MeshValueCollection< T > &mesh)
Definition MeshFunction.h:431
std::shared_ptr< const Mesh > mesh() const
Definition MeshFunction.h:491
std::size_t dim() const
Definition MeshFunction.h:498
MeshFunction(std::shared_ptr< const Mesh > mesh, std::size_t dim, const MeshDomains &domains)
Definition MeshFunction.h:371
void set_all(const T &value)
Definition MeshFunction.h:642
MeshFunction(const MeshFunction< T > &f)
Definition MeshFunction.h:408
bool empty() const
Definition MeshFunction.h:504
~MeshFunction()
Destructor.
Definition MeshFunction.h:126
T * values()
Definition MeshFunction.h:522
void init(std::shared_ptr< const Mesh > mesh, std::size_t dim, std::size_t size)
Definition MeshFunction.h:608
MeshFunction(std::shared_ptr< const Mesh > mesh, const std::string filename)
Definition MeshFunction.h:351
std::size_t size() const
Definition MeshFunction.h:510
std::string str(bool verbose) const
Definition MeshFunction.h:665
MeshFunction(std::shared_ptr< const Mesh > mesh, std::size_t dim, const T &value)
Definition MeshFunction.h:342
void init(std::size_t dim)
Definition MeshFunction.h:572
MeshFunction()
Create empty mesh function.
Definition MeshFunction.h:319
MeshFunction< T > & operator=(const MeshFunction< T > &f)
Definition MeshFunction.h:416
void set_value(std::size_t index, const T &value)
Definition MeshFunction.h:626
T & operator[](const MeshEntity &entity)
Definition MeshFunction.h:528
MeshFunction(std::shared_ptr< const Mesh > mesh)
Definition MeshFunction.h:325
void set_values(const std::vector< T > &values)
Definition MeshFunction.h:634
void init(std::size_t dim, std::size_t size)
Definition MeshFunction.h:586
MeshFunction(std::shared_ptr< const Mesh > mesh, const MeshValueCollection< T > &value_collection)
Definition MeshFunction.h:361
void set_value(std::size_t index, const T &value, const Mesh &mesh)
Compatibility function for use in SubDomains.
Definition MeshFunction.h:261
std::vector< std::size_t > where_equal(T value)
Definition MeshFunction.h:650
void init(std::shared_ptr< const Mesh > mesh, std::size_t dim)
Definition MeshFunction.h:599
const T * values() const
Definition MeshFunction.h:516
MeshFunction(std::shared_ptr< const Mesh > mesh, std::size_t dim)
Definition MeshFunction.h:333
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