librheolef-dev 5.93-2 / usr / include / rheolef / geo-connectivity.h

#ifndef _RHEO_GEO_CONNECTIVITY_H
#define _RHEO_GEO_CONNECTIVITY_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================

// ============================================================================
//  connectivity manipulations
// ============================================================================

#include "rheolef/georep.h"
namespace rheolef { 

// can be used to build explicitely edges:
void build_edge (const geo& omega, std::vector<geo_element>& edge);
void build_edge (const geo& omega, std::vector<std::pair<size_t,size_t> >& edge);

// the rest is used internaly by the class geo and error_estimator_zz

// a := a inter b
// NOTE: could be in O(a.size + b.size) since a and b are sorted
template<class Set>
static
void
in_place_set_intersection (Set& a, const Set& b)
{
    // c := a inter b
    Set c;
    set_intersection (a.begin(), a.end(), b.begin(), b.end(), 
	inserter(c, c.end()));

    // a := c
    a.erase(a.begin(), a.end());
    copy (c.begin(), c.end(), inserter(a, a.end()));
}
// a := a union b
// NOTE: could be in O(a.size + b.size) since a and b are sorted
template<class Set>
static
void
in_place_set_union (Set& a, const Set& b)
{
    // c := a inter b
    Set c;
    set_union (a.begin(), a.end(), b.begin(), b.end(), 
	inserter(c, c.end()));

    // a := c
    a.erase(a.begin(), a.end());
    copy (c.begin(), c.end(), inserter(a, a.end()));
}
template<class Iterator, class SetRandomIterator>
static 
void
build_point_to_element_sets (
  // input
	Iterator           iter,	// elements
	Iterator           last,
  // output
	SetRandomIterator  ball)	// point-to-element-set
{
    typedef typename std::iterator_traits<Iterator>::value_type element_type;
    typedef typename element_type::size_type size_type;

    for (size_type idx = 0; iter != last; iter++, idx++) {
	const element_type& E = *iter;
	for (size_type i = 0; i < E.size(); i++) {
	    ball [E[i]].insert(ball [E[i]].end(), idx);
        }
    }
}
template<class Element, class SetRandomIterator, 
	class Set, class RandomIterator>
void
build_set_that_contains_S (
	const Element&       S,
	SetRandomIterator    ball,
	RandomIterator	     element,
	Set&                 contains_S)
  //! builds a set of elements that all contain S.
{
    typedef geo_element::size_type size_type;

    contains_S.erase(contains_S.begin(), contains_S.end());
    if (S.size() == 0) {
	return;
    }
    copy (ball[S[0]].begin(), ball[S[0]].end(),
	inserter(contains_S, contains_S.end()));

    for (size_type i = 1; i < S.size(); i++) {
        in_place_set_intersection (contains_S, ball[S[i]]);
    }
}
template<class RandomIterator, 
	class RandomIterator2, class RandomIterator3,
	class SetRandomIterator, class Size>
Size
build_subgeo_numbering (
  // input
	RandomIterator           iter,
	RandomIterator           last,
	SetRandomIterator        ball,
	Size			 subgeo_dim,
	RandomIterator2          count_geo,
	RandomIterator3          count_element)
{
    typedef typename std::iterator_traits<RandomIterator>::value_type    element_type;
    typedef reference_element::size_type                   size_type;
    typedef std::set<size_type>				   set_type;
    typedef std::set<size_type>::const_iterator            const_iterator_set;

    RandomIterator element = iter;
    Size n_subgeo = 0;
    for (;iter != last; iter++) {
	const element_type& K = *iter;
	size_type n = K.n_subgeo(subgeo_dim);
	for (size_type j = 0; j < n; j++) {
	    if (K.subgeo(subgeo_dim,j) != geo_element::not_set) {
		// subgeo already numbered in K
		continue;
	    }
	    // found a new subgeo
	    element_type S;
	    K.build_subgeo(subgeo_dim, j, S);
	    count_geo     [S.dimension()] ++;
	    count_element [S.type()] ++;

	    set_type contains_S;
	    build_set_that_contains_S (S, ball, element, contains_S);

	    for (const_iterator_set p = contains_S.begin();
		  p != contains_S.end(); p++) {
		element [*p].set_subgeo(S, size_type(n_subgeo));
	    }
	    n_subgeo++;
	}
    }
    return n_subgeo;
}
template<class RandomIterator, class SetRandomIterator, 
	class RandomIterator2, class RandomIterator3,
	class Size1, class Size2>
std::istream&
load_subgeo_numbering (
  // input
	RandomIterator           element,
	SetRandomIterator        ball,
	Size1			 subgeo_dim,
	Size2			 n_subgeo,
  // in-output	
	std::istream&		 is,
	RandomIterator2          count_geo,
	RandomIterator3          count_element)
{
    typedef typename std::iterator_traits<RandomIterator>::value_type    element_type;
    typedef reference_element::size_type                   size_type;
    typedef std::set<size_type>				   set_type;
    typedef std::set<size_type>::const_iterator            const_iterator_set;

    for (size_type idx = 0; idx < n_subgeo; idx++) {

	element_type S;
	is >> S;
	count_geo     [S.dimension()] ++;
	count_element [S.type()] ++;

	set_type contains_S;
	build_set_that_contains_S (S, ball, element, contains_S);

	for (const_iterator_set p = contains_S.begin();
		  p != contains_S.end(); p++) {
	    element [*p].set_subgeo(S, idx);
	}
    }
    return is;
}
// NOTE: c'est pour les domaines
//  ca pourrait etre fait au vol
//  pendant la lecture
//
template<class RandomIterator, class Iterator,
	class SetRandomIterator>
void
propagate_subgeo_numbering (
  // input
	Iterator          	 iter,
	Iterator           	 last,
	RandomIterator           element,
	SetRandomIterator        ball)
{
    typedef typename std::iterator_traits<Iterator>::value_type    element_type;
    typedef reference_element::size_type                   size_type;
    typedef std::set<size_type>				   set_type;
    typedef std::set<size_type>::const_iterator            const_iterator_set;

    for (; iter != last; iter++) {

	element_type& S = *iter;
	set_type contains_S;
	build_set_that_contains_S (S, ball, element, contains_S);

	if (contains_S.size() == 0) {
	    error_macro ("domain element not in mesh: " << S);
	    return;
	}
	// now, choose any element K containing S
	const_iterator_set p = contains_S.begin();
	element_type& K = element [*p];

	// and propagate numbering from K into S...
        for (size_type d = 1; d <= S.dimension(); d++) {
	    size_type nsubgeo = S.n_subgeo(d);
	    for (size_type i = 0; i < nsubgeo; i++) {
		element_type E;
                S.build_subgeo(d, i, E);
		size_type idx = K.subgeo(E);
		S.set_subgeo(d, i, idx);
	    }
	}
    }
}
template<class RandomIterator, class SetRandomIterator,
	class Container, class Size>
void
build_boundary (
  // input
	RandomIterator           iter,
	RandomIterator           last,
	SetRandomIterator        ball,
	Size			 d,
  // output:
        Container&		 boundary)
{
    typedef typename std::iterator_traits<RandomIterator>::value_type    element_type;
    typedef reference_element::size_type                   size_type;
    typedef std::set<size_type>				   set_type;
    typedef std::set<size_type>::const_iterator            const_iterator_set;

    RandomIterator element = iter;
    for (;iter != last; iter++) {
	const element_type& K = *iter;
	size_type n = K.n_subgeo(d);
	for (size_type j = 0; j < n; j++) {
	    element_type S;
	    K.build_subgeo(d, j, S);
	    S.set_index(K.subgeo(d,j));
	    set_type contains_S;
	    build_set_that_contains_S (S, ball, element, contains_S);
	    switch (contains_S.size()) {
		case 1:
		    boundary.push_back (S);
		    break;
		case 2:
		    // interior face
		    break;
		default:
		    error_macro ("face " << S << " shared by " << contains_S.size() << " elements");
	    }
	}
    }
}
void build_connectivity (georep& g);

}// namespace rheolef
#endif // _RHEO_GEO_CONNECTIVITY_H