level_set.cc
Go to the documentation of this file.
#include "rheolef/level_set.h"
namespace rheolef {
// --------------------------------------------------------------------------------
// gestion de la numerotation locale
// --------------------------------------------------------------------------------
// TODO: move in reference_element
static
size_t
edge_t_iloc (size_t l, size_t m)
{
static const int edge_t_iloc_table [3][3] = {
{-1, 0, 2},
{ 0,-1, 1},
{ 2, 1,-1}};
return size_t(edge_t_iloc_table [l][m]);
}
static
size_t
edge_T_iloc (size_t l, size_t m)
{
static const int edge_T_iloc_table [4][4] = {
{-1, 0, 2, 3},
{ 0,-1, 1, 4},
{ 2, 1,-1, 5},
{ 3, 4, 5,-1}};
return size_t(edge_T_iloc_table [l][m]);
}
static
size_t
face_T_iloc (size_t l, size_t m, size_t n)
{
static size_t face_T_iloc_table [4][4][4];
bool static initialized = false;
if (!initialized) {
for (size_t i = 0; i < 4; i++)
for (size_t j = 0; j < 4; j++)
for (size_t k = 0; k < 4; k++)
face_T_iloc_table [i][j][k] = size_t(-1);
reference_element hat_K (reference_element::T);
for (size_t i_face = 0; i_face < hat_K.n_face(); i_face++) {
size_t p[3];
for (size_t k = 0; k < 3; k++) p[k] = hat_K.subgeo_local_vertex(2,i_face,k);
face_T_iloc_table [p[0]][p[1]][p[2]] = i_face;
face_T_iloc_table [p[0]][p[2]][p[1]] = i_face;
face_T_iloc_table [p[1]][p[0]][p[2]] = i_face;
face_T_iloc_table [p[1]][p[2]][p[0]] = i_face;
face_T_iloc_table [p[2]][p[0]][p[1]] = i_face;
face_T_iloc_table [p[2]][p[1]][p[0]] = i_face;
}
}
return face_T_iloc_table [l][m][n];
}
// --------------------------------------------------------------------------------
// gestion de la precision
// --------------------------------------------------------------------------------
static Float level_set_epsilon = 100*std::numeric_limits<Float>::epsilon();
template <class T>
static
bool
is_zero (const T& x) {
// TODO: control by level_set_option ?
return fabs(x) <= level_set_epsilon;
}
template <class T>
static
bool
have_same_sign (const T& x, const T& y) {
using namespace details;
return !is_zero(x) && !is_zero(y) && x*y > 0;
}
template <class T>
static
bool
have_opposite_sign (const T& x, const T& y) {
using namespace details;
return !is_zero(x) && !is_zero(y) && x*y < 0;
}
// --------------------------------------------------------------------------------
// 2D: fonctions locales : sur un seul element triangle
// --------------------------------------------------------------------------------
// appele lors du 1er passage qui liste les elements de la bande cas de dimension 2
template <class T>
static
bool
belongs_to_band_t (const std::vector<T>& f) {
using namespace details;
if (have_same_sign(f[0],f[1]) && have_same_sign(f[0],f[2])) return false;
// on rejette le triangle dans tous les sommets de meme signe :
if (is_zero(f[0]) && is_zero(f[1]) && is_zero(f[2])) return false;
// on rejette les triangles dont un sommet :
if (is_zero(f[0]) && have_same_sign(f[1],f[2])) return false;
if (is_zero(f[1]) && have_same_sign(f[0],f[2])) return false;
if (is_zero(f[2]) && have_same_sign(f[0],f[1])) return false;
return true;
}
// apellee lors du calcul des matrices M_K et A_K pour K dans la bande cas de dimension 2
template <class T>
static
size_t
isolated_vertex_t (const std::vector<T>& f)
{
using namespace details;
/* on retourne le sommet isole a chaque fois */
if (have_same_sign (f[0],f[1]) && have_opposite_sign(f[0],f[2])) return 2;
if (have_opposite_sign(f[0],f[1]) && have_same_sign (f[0],f[2])) return 1;
if (have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2])) return 0;
if (is_zero(f[0]) && have_opposite_sign(f[1],f[2])) return 1; /* on peut retourner 2 de meme*/
if (is_zero(f[1]) && have_opposite_sign(f[0],f[2])) return 0; /* on peut retourner 2 */
if (is_zero(f[2]) && have_opposite_sign(f[0],f[1])) return 0; /* on peut retourner 1 */
if (is_zero(f[0]) && is_zero(f[1]) && ! is_zero(f[2])) return 2;
if (is_zero(f[0]) && is_zero(f[2]) && ! is_zero(f[1])) return 1;
return 0; /* f1 == 0 et f2 == 0 et f0 != 0 */
}
template <class T>
static
void
subcompute_matrix_t (
const std::vector<point_basic<T> >& x,
const std::vector<T>& f,
std::vector<size_t>& j,
point_basic<T>& a,
point_basic<T>& b,
T& S)
{
using namespace details;
j.resize (3);
j[0] = isolated_vertex_t (f);
j[1] = (j[0]+1) % 3;
j[2] = (j[1]+1) % 3;
// edge {j1,j2} has normal oriented as grad(f), in f>0 direction:
if (! is_zero(f[j[0]]) && f[j[0]] < 0) std::swap (j[1], j[2]);
T theta_1= f[j[1]]/(f[j[1]]-f[j[0]]);
T theta_2= f[j[2]]/(f[j[2]]-f[j[0]]);
// calcul des coordonnes d'intersection
a = theta_1*x[j[0]]+(1-theta_1)*x[j[1]];
b = theta_2*x[j[0]]+(1-theta_2)*x[j[2]];
S = sqrt(pow(a[0]-b[0],2)+pow(a[1]-b[1],2));
if (is_zero(f[j[1]]) && is_zero(f[j[2]])) {
S /= 2;
}
}
// --------------------------------------------------------------------------------
// 3D: fonctions locales : sur un seul element tetraedre
// --------------------------------------------------------------------------------
class quadruplet {
public:
quadruplet (size_t a=0, size_t b=0, size_t c=0, size_t d=0) {
q[0]=a;
q[1]=b;
q[2]=c;
q[3]=d;
}
size_t operator[] (size_t i) const {
return q[i%4];
}
size_t& operator[] (size_t i) {
return q[i%4];
}
friend std::ostream& operator<< (std::ostream& out, const quadruplet& q) {
out << "((" << q[0] << "," << q[1] << "), (" << q[2] << "," << q[3] << "))";
return out;
}
protected:
size_t q[4];
};
// appele lors du 1er passage qui liste les elements de la bande cas de dimension 3
template <class T>
static
bool
belongs_to_band_T (const std::vector<T>& f)
{
using namespace details;
if (have_same_sign(f[0],f[1]) && have_same_sign(f[0],f[2]) && have_same_sign(f[2],f[3])) return false;
// cas ou 4 points s'annulent en dimension 3 est degenere
if (is_zero(f[0]) && is_zero(f[1]) && is_zero(f[2]) && is_zero(f[3])) return false;
if (is_zero(f[0]) && have_same_sign(f[1],f[2]) && have_same_sign(f[1],f[3])) return false;
if (is_zero(f[1]) && have_same_sign(f[0],f[2]) && have_same_sign(f[0],f[3])) return false;
if (is_zero(f[2]) && have_same_sign(f[0],f[1]) && have_same_sign(f[0],f[3])) return false;
if (is_zero(f[3]) && have_same_sign(f[0],f[1]) && have_same_sign(f[0],f[2])) return false;
// cas ou f s'annule sur 2 sommets et garde le meme signe sur les 2 autres sommets est exclu
if (is_zero(f[0]) && is_zero(f[1]) && have_same_sign(f[2],f[3])) return false;
if (is_zero(f[0]) && is_zero(f[2]) && have_same_sign(f[1],f[3])) return false;
if (is_zero(f[0]) && is_zero(f[3]) && have_same_sign(f[1],f[2])) return false;
if (is_zero(f[1]) && is_zero(f[2]) && have_same_sign(f[0],f[3])) return false;
if (is_zero(f[1]) && is_zero(f[3]) && have_same_sign(f[0],f[2])) return false;
if (is_zero(f[2]) && is_zero(f[3]) && have_same_sign(f[1],f[0])) return false;
return true;
}
// apellee lors du calcul des matrices M_K et A_K pour T dans la bande cas de dimension
template <class T>
bool
intersection_is_quadrilateral_T (const std::vector<T>& f, quadruplet& q)
{
if (have_same_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2]) && have_same_sign(f[2],f[3])) {
if (f[0] > 0) q = quadruplet(0,1, 2,3);
else q = quadruplet(2,3, 0,1);
return true;
}
if (have_opposite_sign(f[0],f[1]) && have_same_sign(f[0],f[2]) && have_opposite_sign(f[2],f[3])) {
if (f[0] < 0) q = quadruplet(0,2, 1,3);
else q = quadruplet(1,3, 0,2);
return true;
}
if (have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2]) && have_opposite_sign(f[2],f[3])) {
if (f[0] > 0) q = quadruplet(0,3, 1,2);
else q = quadruplet(1,2, 0,3);
return true;
}
return false;
}
// cas d'une intersection triangle:
template <class T>
static
size_t
isolated_vertex_T (const std::vector<T>& f)
{
using namespace details;
// cas ou l'intersection est un triangle
if (have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2]) && have_same_sign (f[2],f[3])) return 0;
if (have_same_sign (f[0],f[1]) && have_opposite_sign(f[0],f[2]) && have_opposite_sign(f[2],f[3])) return 2;
if (have_same_sign (f[0],f[1]) && have_same_sign (f[0],f[2]) && have_opposite_sign(f[2],f[3])) return 3;
// cas ou f s'annule sur un sommet et change de signe sur les 2 autres sommets
if (have_opposite_sign(f[0],f[1]) && have_same_sign(f[0],f[2]) && have_same_sign(f[2],f[3])) return 1;
if (is_zero(f[0]) && have_same_sign (f[1],f[2]) && have_opposite_sign(f[1],f[3])) return 3;
if (is_zero(f[0]) && have_opposite_sign(f[1],f[2]) && have_same_sign (f[1],f[3])) return 2;
if (is_zero(f[0]) && have_opposite_sign(f[1],f[2]) && have_opposite_sign(f[1],f[3])) return 1;
if (is_zero(f[1]) && have_opposite_sign(f[0],f[2]) && have_same_sign (f[0],f[3])) return 2;
if (is_zero(f[1]) && have_same_sign (f[0],f[2]) && have_opposite_sign(f[0],f[3])) return 3;
if (is_zero(f[1]) && have_opposite_sign(f[0],f[2]) && have_opposite_sign(f[0],f[3])) return 0;
if (is_zero(f[2]) && have_opposite_sign(f[0],f[1]) && have_same_sign (f[0],f[3])) return 1;
if (is_zero(f[2]) && have_same_sign (f[0],f[1]) && have_opposite_sign(f[0],f[3])) return 3;
if (is_zero(f[2]) && have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[3])) return 0;
if (is_zero(f[3]) && have_opposite_sign(f[0],f[1]) && have_same_sign (f[0],f[2])) return 1;
if (is_zero(f[3]) && have_same_sign (f[0],f[1]) && have_opposite_sign(f[0],f[2])) return 2;
if (is_zero(f[3]) && have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2])) return 0;
// cas ou f s'annule en 2 sommets et change de signe sur les 2 autres sommets
if (is_zero(f[0]) && is_zero(f[1]) && have_opposite_sign(f[2],f[3])) return 2; // ou 3
if (is_zero(f[0]) && is_zero(f[2]) && have_opposite_sign(f[1],f[3])) return 1;
if (is_zero(f[0]) && is_zero(f[3]) && have_opposite_sign(f[1],f[2])) return 1;
if (is_zero(f[1]) && is_zero(f[2]) && have_opposite_sign(f[0],f[3])) return 0;
if (is_zero(f[1]) && is_zero(f[3]) && have_opposite_sign(f[0],f[2])) return 0;
if (is_zero(f[2]) && is_zero(f[3]) && have_opposite_sign(f[0],f[1])) return 0;
// le triangle d'intersection est la face du tetradre ou f s'annule sur les 3 sommets
if (is_zero(f[0]) && is_zero(f[1]) && is_zero(f[2]) && ! is_zero(f[3])) return 3;
if (is_zero(f[0]) && is_zero(f[1]) && is_zero(f[3]) && ! is_zero(f[2])) return 2;
if (is_zero(f[1]) && is_zero(f[2]) && is_zero(f[3]) && ! is_zero(f[0])) return 0;
return 1;
}
template <class T>
static
void
subcompute_matrix_triangle_T (
const std::vector<point_basic<T> >& x,
const std::vector<T>& f,
std::vector<size_t>& j,
point_basic<T>& a,
point_basic<T>& b,
point_basic<T>& c,
T& aire)
{
using namespace details;
j.resize(4);
j[0] = isolated_vertex_T (f);
j[1] = (j[0]+1) % 4;
j[2] = (j[1]+1) % 4;
j[3] = (j[2]+1) % 4;
// orient
if (! is_zero(f[j[0]]) && ((f[j[0]] > 0 && j[0] % 2 == 0) || (f[j[0]] < 0 && j[0] % 2 == 1)))
std::swap (j[1], j[2]);
T theta_1 = f[j[1]]/(f[j[1]]-f[j[0]]);
T theta_2 = f[j[2]]/(f[j[2]]-f[j[0]]);
T theta_3 = f[j[3]]/(f[j[3]]-f[j[0]]);
/* calcul des coordonnees d'intersection */
a = theta_1*x[j[0]]+(1-theta_1)*x[j[1]];
b = theta_2*x[j[0]]+(1-theta_2)*x[j[2]];
c = theta_3*x[j[0]]+(1-theta_3)*x[j[3]];
aire = 0.5* norm (vect( b-a, c-a));
if (is_zero(f[j[1]]) && is_zero(f[j[2]]) && is_zero(f[j[3]])) {
aire /= 2;
}
}
template <class T>
static
void
subcompute_matrix_quadrilateral_T (
const std::vector<point_basic<T> >& x,
const std::vector<T>& f,
const quadruplet& q,
point_basic<T>& a,
point_basic<T>& b,
point_basic<T>& c,
point_basic<T>& d,
T& aire_Q)
{
// intersection:
// a = segment {x(q0) x(q2)} inter {f=0}
// b = segment {x(q1) x(q2)} inter {f=0}
// d = segment {x(q1) x(q3)} inter {f=0}
// c = segment {x(q0) x(q3)} inter {f=0}
// quadrilatere abdc = triangle(abd) union triangle(adc)
T theta_1 = f[q[2]]/(f[q[2]]-f[q[0]]);
T theta_2 = f[q[2]]/(f[q[2]]-f[q[1]]);
T theta_3 = f[q[3]]/(f[q[3]]-f[q[0]]);
T theta_4 = f[q[3]]/(f[q[3]]-f[q[1]]);
/* calcul des coordonnees d'intersection */
a = theta_1*x[q[0]]+(1-theta_1)*x[q[2]];
b = theta_2*x[q[1]]+(1-theta_2)*x[q[2]];
c = theta_3*x[q[0]]+(1-theta_3)*x[q[3]];
d = theta_4*x[q[1]]+(1-theta_4)*x[q[3]];
aire_Q = 0.5*norm(vect(a-c,a-b)) + 0.5*norm(vect(d-c,d-b));
}
// --------------------------------------------------------------------------------
// level_set: compte the intersection mesh
// --------------------------------------------------------------------------------
typedef geo_element_auto<heap_allocator<geo_element::size_type> > element_type;
template <class T, class M>
void
gamma_list2disarray (
const std::list<point_basic<T> >& gamma_node_list,
std::array<std::list<std::pair<element_type,size_t> >,
reference_element::max_variant> gamma_side_list,
const communicator& comm,
size_t d,
disarray<point_basic<T>, M>& gamma_node,
std::array<disarray<element_type,M>,
reference_element::max_variant>& gamma_side,
disarray<size_t,M>& sid_ie2bnd_ie)
{
typedef geo_element::size_type size_type;
// 1) nodes:
size_type nnod = gamma_node_list.size();
distributor gamma_node_ownership (distributor::decide, comm, nnod);
gamma_node.resize (gamma_node_ownership);
typename disarray<point_basic<T>, M>::iterator node_iter = gamma_node.begin();
for (typename std::list<point_basic<T> >::const_iterator
iter = gamma_node_list.begin(),
last = gamma_node_list.end();
iter != last; iter++, node_iter++) {
*node_iter = *iter;
}
// 2) sides:
heap_allocator<size_type> alloc;
size_type map_dim = d-1;
size_type order = 1;
size_type nsid = 0;
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
size_type nsidv = gamma_side_list [variant].size();
distributor gamma_sidv_ownership (distributor::decide, comm, nsidv);
nsid += nsidv;
element_type element_init (variant, order, alloc);
gamma_side[variant].resize (gamma_sidv_ownership, element_init);
typename disarray<element_type, M>::iterator side_iter = gamma_side[variant].begin();
for (typename std::list<std::pair<element_type,size_type> >::const_iterator
iter = gamma_side_list[variant].begin(),
last = gamma_side_list[variant].end();
iter != last; iter++, side_iter++) {
*side_iter = (*iter).first;
}
}
// 3) side2band correspondance
distributor gamma_sid_ownership (distributor::decide, comm, nsid);
const size_type undef = std::numeric_limits<size_t>::max();
sid_ie2bnd_ie.resize (gamma_sid_ownership, undef);
typename disarray<size_type, M>::iterator idx_iter = sid_ie2bnd_ie.begin();
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
for (typename std::list<std::pair<element_type,size_type> >::const_iterator
iter = gamma_side_list[variant].begin(),
last = gamma_side_list[variant].end();
iter != last; iter++, idx_iter++) {
*idx_iter = (*iter).second;
}
}
}
struct to_solve {
size_t variant, S_ige, k, dis_i;
to_solve (size_t variant1, size_t S_ige1=0, size_t k1=0, size_t dis_i1=0)
: variant(variant1), S_ige(S_ige1), k(k1), dis_i(dis_i1) {}
};
template <class T, class M>
geo_basic<T,M>
level_set_internal (
const field_basic<T,M>& fh,
const level_set_option& opt,
std::vector<size_t>& bnd_dom_ie_list,
disarray<size_t,M>& sid_ie2bnd_ie)
{
using namespace std;
using namespace details;
typedef geo_element::size_type size_type;
level_set_epsilon = opt.epsilon; // set global variable
fh.dis_dof_update();
const geo_basic<T,M>& lambda = fh.get_geo();
const space_basic<T,M>& Xh = fh.get_space();
check_macro(lambda.order() == 1, "Only order=1 level set mesh supported");
check_macro(fh.get_approx() == "P1", "Only P1 level set function supported");
size_type order = 1;
std::vector<size_type> dis_idof;
std::vector<T> f;
size_type d = lambda.dimension();
size_type map_dim = d-1;
heap_allocator<size_type> alloc;
std::array<list<pair<element_type,size_type> >,
reference_element::max_variant> gamma_side_list;
list<point_basic<T> > gamma_node_list;
std::vector<size_type> j(d+1);
std::vector<point_basic<T> > x(d+1);
const size_type not_marked = numeric_limits<size_t>::max();
const size_type undef = numeric_limits<size_t>::max();
distributor node_ownership = lambda.sizes().node_ownership;
size_type first_dis_inod = node_ownership.first_index();
disarray<size_type> marked_node (node_ownership, not_marked);
disarray<size_type> extern_node (node_ownership, not_marked);
set<size_type> ext_marked_node_idx;
list<to_solve> node_to_solve;
distributor edge_ownership = lambda.sizes().ownership_by_dimension[1];
size_type first_dis_iedg = edge_ownership.first_index();
disarray<size_type> marked_edge (edge_ownership, not_marked);
disarray<std::pair<size_type,point_basic<T> > >
extern_edge (edge_ownership, std::make_pair(not_marked,point_basic<T>()));
set<size_type> ext_marked_edge_idx;
list<to_solve> edge_to_solve;
distributor face_ownership = lambda.sizes().ownership_by_dimension[2];
size_type first_dis_ifac = face_ownership.first_index();
disarray<size_type> marked_face (face_ownership, not_marked);
communicator comm = node_ownership.comm();
size_type my_proc = comm.rank();
bnd_dom_ie_list.resize(0);
// ------------------------------------------------------------
// 1) loop on lambda & build intersection sides
// ------------------------------------------------------------
for (size_type ie = 0, ne = lambda.size(), bnd_ie = 0; ie < ne; ie++) {
// ---------------------------------------------------------
// 1.1) fast check if there is an intersection:
// ---------------------------------------------------------
const geo_element& K = lambda [ie];
Xh.dis_idof (K, dis_idof);
f.resize (dis_idof.size());
for (size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
f [loc_idof] = fh.dis_dof (dis_idof[loc_idof]);
}
bool do_intersect = false;
switch (K.variant()) {
case reference_element::t:
do_intersect = belongs_to_band_t (f);
break;
case reference_element::T: {
do_intersect = belongs_to_band_T (f);
break;
}
default :
error_macro("level set: element type `" << K.name() << "' not yet supported");
}
if (! do_intersect) continue;
bnd_dom_ie_list.push_back (ie);
// ---------------------------------------------------------
// 1.2) compute the intersection
// ---------------------------------------------------------
x.resize (dis_idof.size());
for (size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
size_type loc_inod = loc_idof; // assume here that fh has an isoparametric approx
size_type dis_inod = K [loc_inod];
x [loc_idof] = lambda.dis_node(dis_inod);
}
element_type S (alloc);
switch (K.variant()) {
case reference_element::t: {
// ---------------------------------------------------------
// 1.2.a) triangle -> 1 edge
// ---------------------------------------------------------
point_basic<T> a, b;
T length;
subcompute_matrix_t (x, f, j, a, b, length);
if (is_zero(f[j[1]]) && is_zero(f[j[2]])) {
// ---------------------------------------------------------
// 1.2.a.1) edge {j1,j2} is included in the bbox mesh
// ---------------------------------------------------------
for (size_type k = 0; k < 2; k++) {
size_type dis_inod = K[j[k+1]];
if (node_ownership.is_owned (dis_inod)) {
size_type inod = dis_inod - first_dis_inod;
if (marked_node [inod] == not_marked) {
marked_node [inod] = gamma_node_list.size();
gamma_node_list.push_back (lambda.node(inod));
}
} else {
// dis_inod is owned by another partition jproc:
// if there is another neighbour element K' from the same partition jproc
// then jproc will insert dis_inod in gamma
// so, there is nothing to do here
// otherwise, dis_inod will be orphan in jproc and will be re-affected to my_proc
extern_node.dis_entry (dis_inod) = my_proc;
}
}
size_type loc_iedg = edge_t_iloc (j[1], j[2]);
size_type dis_iedg = K.edge (loc_iedg);
if (edge_ownership.is_owned (dis_iedg)) {
size_type iedg = dis_iedg - first_dis_iedg;
if (marked_edge [iedg] == not_marked) {
S.reset(reference_element::e, order);
marked_edge [iedg] = gamma_side_list [S.variant()].size();
size_type S_ige = marked_edge [iedg];
for (size_type k = 0; k < S.n_node(); k++) {
size_type dis_inod = K[j[k+1]];
if (node_ownership.is_owned (dis_inod)) {
size_type inod = dis_inod - first_dis_inod;
S[k] = marked_node [inod];
} else {
S[k] = undef;
node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
ext_marked_node_idx.insert (dis_inod);
}
}
gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
}
} else {
// intersection is an edge of lambda that is owned by another partition jproc
// the neighbour element K' across this e"dge belongs to the same partition jproc
// and will insert dis_iedg in gamma
// so, there is nothing to do here
}
} else {
// ---------------------------------------------------------
// 1.2.a.2) edge {j1,j2} is interior to the triangle
// ---------------------------------------------------------
S.reset(reference_element::e, order);
size_type S_ige = gamma_side_list [S.variant()].size ();
point_basic<T> xx[2] = {a,b};
for (size_type k = 0; k < 2; k++) {
if (! is_zero(f[j[k+1]]) && ! is_zero(f[j[0]])) {
// xk is inside edge {j0,j[k+1]} of triangle K:
size_type loc_iedg = edge_t_iloc (j[0], j[k+1]);
size_type dis_iedg = K.edge (loc_iedg);
if (edge_ownership.is_owned (dis_iedg)) {
size_type iedg = dis_iedg - first_dis_iedg;
if (marked_edge [iedg] == not_marked) {
marked_edge [iedg] = gamma_node_list.size();
gamma_node_list.push_back (xx[k]);
}
S[k] = marked_edge [iedg];
} else {
S[k] = undef;
edge_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_iedg));
ext_marked_edge_idx.insert (dis_iedg);
extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
}
} else { // xk is at edge boundary: a node of the 2d mesh
size_type dis_inod = (!is_zero(f[j[0]])) ? K[j[k+1]] : K[j[0]];
if (node_ownership.is_owned (dis_inod)) {
size_type inod = dis_inod - first_dis_inod;
if (marked_node [inod] == not_marked) {
marked_node [inod] = gamma_node_list.size();
gamma_node_list.push_back (lambda.node(inod));
}
S[k] = marked_node [inod];
} else {
S[k] = undef;
node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
ext_marked_node_idx.insert (dis_inod);
extern_node.dis_entry (dis_inod) = my_proc;
}
}
}
// S[0] == S[1] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
check_macro (S[0] != S[1] || S[0] == undef, "degenerate 2d intersection");
gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
}
break;
}
case reference_element::T: {
// ---------------------------------------------------------
// 1.2.b) tetrahedron -> 1 triangle or 1 quadrangle
// ---------------------------------------------------------
quadruplet q;
point_basic<T> a, b, c, d;
T aire;
if (!intersection_is_quadrilateral_T (f, q)) {
subcompute_matrix_triangle_T (x, f, j, a, b, c, aire);
if (is_zero(f[j[1]]) && is_zero(f[j[2]]) && is_zero(f[j[3]])) {
// the full face {j1,j2,j3} is included in the surface mesh:
for (size_type k = 0; k < 3; k++) {
size_type dis_inod = K[j[k+1]];
if (node_ownership.is_owned (dis_inod)) {
size_type inod = dis_inod - first_dis_inod;
if (marked_node [inod] == not_marked) {
marked_node [inod] = gamma_node_list.size();
gamma_node_list.push_back (lambda.node(inod));
}
} else {
// dis_inod is owned by another partition jproc:
// if there is another neighbour element K' from the same partition jproc
// then jproc will insert dis_inod in gamma
// so, there is nothing to do here
// otherwise, dis_inod will be orphan in jproc and will be re-affected to my_proc
extern_node.dis_entry (dis_inod) = my_proc;
}
}
size_type loc_ifac = face_T_iloc (j[1], j[2], j[3]);
size_type dis_ifac = K.face (loc_ifac);
if (face_ownership.is_owned (dis_ifac)) {
size_type ifac = dis_ifac - first_dis_ifac;
if (marked_face [ifac] == not_marked) {
S.reset(reference_element::t, order);
marked_face [ifac] = gamma_side_list [S.variant()].size();
size_type S_ige = marked_face [ifac];
for (size_type k = 0; k < S.n_node(); k++) {
size_type dis_inod = K[j[k+1]];
if (node_ownership.is_owned (dis_inod)) {
size_type inod = dis_inod - first_dis_inod;
S[k] = marked_node [inod];
} else {
S[k] = undef;
node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
ext_marked_node_idx.insert (dis_inod);
}
}
gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
} else {
// the side will be inserted by the neighbour of K in another partition
// so, nothing to do
}
} else {
// intersection is a face of lambda that is owned by another partition jproc
// the neighbour element K' across this face belongs to the same partition jproc
// and will insert dis_ifac in gamma
// so, there is nothing to do here
}
} else {
// create the new face {j1,j2,j3} by intersections:
S.reset(reference_element::t, order);
size_type S_ige = gamma_side_list [S.variant()].size();
point_basic<T> xx[3] = {a,b,c};
for (size_type k = 0; k < 3; k++) {
if (! is_zero(f[j[k+1]]) && ! is_zero(f[j[0]])) {
// xk is inside edge {j0,j[k+1]} of triangle K:
size_type loc_iedg = edge_T_iloc (j[0], j[k+1]);
size_type dis_iedg = K.edge (loc_iedg);
if (edge_ownership.is_owned (dis_iedg)) {
size_type iedg = dis_iedg - first_dis_iedg;
if (marked_edge [iedg] == not_marked) {
marked_edge [iedg] = gamma_node_list.size();
gamma_node_list.push_back (xx[k]);
}
S[k] = marked_edge [iedg];
} else {
S[k] = undef;
edge_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_iedg));
ext_marked_edge_idx.insert (dis_iedg);
extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
}
} else { // xk is at edge boundary: a node of the 2d mesh
size_type dis_inod = (!is_zero(f[j[0]])) ? K[j[k+1]] : K[j[0]];
if (node_ownership.is_owned (dis_inod)) {
size_type inod = dis_inod - first_dis_inod;
if (marked_node [inod] == not_marked) {
marked_node [inod] = gamma_node_list.size();
gamma_node_list.push_back (lambda.node(inod));
}
S[k] = marked_node [inod];
} else {
S[k] = undef;
node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
ext_marked_node_idx.insert (dis_inod);
extern_node.dis_entry (dis_inod) = my_proc;
}
}
}
// S[0] == S[j] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
check_macro ((S[0] != S[1] || S[0] == undef) &&
(S[1] != S[2] || S[1] == undef) &&
(S[2] != S[0] || S[2] == undef), "degenerate 3d intersection");
gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
}
} else {
// create the new quadri face by intersections:
subcompute_matrix_quadrilateral_T (x, f, q, a, b, c, d, aire);
S.reset (reference_element::q, order);
size_type S_ige = gamma_side_list[reference_element::q].size();
size_type S1_ige = gamma_side_list[reference_element::t].size(); // or S1+S2=two 't'
size_type S2_ige = S1_ige + 1;
size_type iloc_S1[4] = {0, 1, 2, undef}; // the way to slip a q into 2 t
size_type iloc_S2[4] = {0, undef, 1, 2};
point_basic<T> xx[4] = {a,b,d,c};
size_type s[4] = {q[0],q[2],q[1],q[3]};
for (size_type k = 0; k < 4; k++) {
size_type k1 = (k+1) % 4;
if (! is_zero(f[s[k]]) && ! is_zero(f[s[k1]])) {
// xk is inside edge {j0,j[k+1]} of triangle K:
size_type loc_iedg = edge_T_iloc (s[k], s[k1]);
size_type dis_iedg = K.edge (loc_iedg);
if (edge_ownership.is_owned (dis_iedg)) {
size_type iedg = dis_iedg - first_dis_iedg;
if (marked_edge [iedg] == not_marked) {
marked_edge [iedg] = gamma_node_list.size();
gamma_node_list.push_back (xx[k]);
}
S[k] = marked_edge [iedg];
} else {
S[k] = undef;
if (!opt.split_to_triangle) {
edge_to_solve.push_back (to_solve (reference_element::q, S_ige, k, dis_iedg));
} else {
if (iloc_S1[k] != undef) edge_to_solve.push_back (to_solve (reference_element::t, S1_ige, iloc_S1[k], dis_iedg));
if (iloc_S2[k] != undef) edge_to_solve.push_back (to_solve (reference_element::t, S2_ige, iloc_S2[k], dis_iedg));
}
ext_marked_edge_idx.insert (dis_iedg);
extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
}
} else { // xk is at edge boundary: a node of the 2d mesh
size_type dis_inod = is_zero(f[s[k]]) ? K[s[k]] : K[s[k1]];
if (node_ownership.is_owned (dis_inod)) {
size_type inod = dis_inod - first_dis_inod;
if (marked_node [inod] == not_marked) {
marked_node [inod] = gamma_node_list.size();
gamma_node_list.push_back (lambda.node(inod));
}
S[k] = marked_node [inod];
} else {
S[k] = undef;
if (!opt.split_to_triangle) {
node_to_solve.push_back (to_solve (reference_element::q, S_ige, k, dis_inod));
} else {
if (iloc_S1[k] != undef) node_to_solve.push_back (to_solve (reference_element::t, S1_ige, iloc_S1[k], dis_inod));
if (iloc_S2[k] != undef) node_to_solve.push_back (to_solve (reference_element::t, S2_ige, iloc_S2[k], dis_inod));
}
ext_marked_node_idx.insert (dis_inod);
}
}
}
if (!opt.split_to_triangle) {
check_macro ((S[0] != S[1] || S[0] == undef) &&
(S[1] != S[2] || S[1] == undef) &&
(S[2] != S[3] || S[2] == undef) &&
(S[3] != S[0] || S[3] == undef), "degenerate 3d intersection");
// S[0] == S[j] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
} else {
// split quadri into 2 triangles
// one K -> two (S1,S2) faces: table element2face may return a pair of size_t
// but S1-> and S2->K only is required during assembly
element_type S1 (alloc);
S1.reset (reference_element::t, order);
for (size_type k = 0; k < 4; k++) {
if (iloc_S1[k] != undef) S1 [iloc_S1[k]] = S[k];
}
check_macro ((S1[0] != S1[1] || S1[0] == undef) &&
(S1[1] != S1[2] || S1[1] == undef) &&
(S1[2] != S1[0] || S1[2] == undef), "degenerate 3d intersection");
gamma_side_list [S1.variant()].push_back (make_pair(S1,bnd_ie));
element_type S2 (alloc);
S2.reset (reference_element::t, order);
for (size_type k = 0; k < 4; k++) {
if (iloc_S2[k] != undef) S2 [iloc_S2[k]] = S[k];
}
check_macro ((S2[0] != S2[1] || S2[0] == undef) &&
(S2[1] != S2[2] || S2[1] == undef) &&
(S2[2] != S2[0] || S2[2] == undef), "degenerate 3d intersection");
gamma_side_list [S2.variant()].push_back (make_pair(S2,bnd_ie));
}
} // if-else
break;
}
default : {
error_macro("level set intersection: element not yet implemented: " << K.name());
}
}
bnd_ie++;
}
extern_node.dis_entry_assembly();
extern_edge.dis_entry_assembly();
// ------------------------------------------------------------
// 2) solve orphan nodes, if any
// ------------------------------------------------------------
// 2.1.a) re-affect orphan node to another process where gamma use it
distributor comm_ownership (comm.size(), comm, 1);
disarray<index_set,M> orphan_node (comm_ownership, index_set());
for (size_type inod = 0, nnod = node_ownership.size(); inod < nnod; inod++) {
if (!(extern_node [inod] != not_marked && marked_node [inod] == not_marked)) continue;
// inod is orphan in this proc: not used for gamma (but used for lambda)
// re-affect it to a process that use it for gamma
size_type iproc = extern_node[inod];
size_type dis_inod = first_dis_inod + inod;
index_set dis_inod_set;
dis_inod_set.insert (dis_inod);
orphan_node.dis_entry (iproc) += dis_inod_set;
}
orphan_node.dis_entry_assembly();
// 2.1.b) re-affect orphan edge to another process where gamma use it
// there could be orphan edge in 2d if lambda is a part of a regular mesh:
// => a boundary edge can belong to another proc. This is not yet handled
disarray<index_set,M> orphan_edge (comm_ownership, index_set());
for (size_type iedg = 0, nedg = edge_ownership.size(); iedg < nedg; iedg++) {
if (!(extern_edge [iedg].first != not_marked && marked_edge [iedg] == not_marked)) continue;
// iedg is orphan in this proc: not used for gamma (but used for lambda)
// re-affect it to a process that use it for gamma
size_type iproc = extern_edge[iedg].first;
size_type dis_iedg = first_dis_iedg + iedg;
index_set dis_iedg_set;
dis_iedg_set.insert (dis_iedg);
orphan_edge.dis_entry (iproc) += dis_iedg_set;
}
orphan_edge.dis_entry_assembly ();
check_macro (orphan_edge[0].size() == 0, "unexpected orphan edges");
// 2.2) count total nodes used by gamma
size_type orphan_gamma_nnod = orphan_node[0].size();
size_type regular_gamma_nnod = gamma_node_list.size();
size_type gamma_nnod = regular_gamma_nnod + orphan_gamma_nnod;
distributor gamma_node_ownership (distributor::decide, comm, gamma_nnod);
// 2.3) shift marked_node & marked_edge from gamma_inod local count to gamma_dis_inod one
size_type gamma_first_dis_inod = gamma_node_ownership.first_index();
for (size_type inod = 0, nnod = node_ownership.size(); inod < nnod; inod++) {
if (marked_node [inod] == not_marked) continue;
marked_node [inod] += gamma_first_dis_inod;
}
for (size_type iedg = 0, nedg = edge_ownership.size(); iedg < nedg; iedg++) {
if (marked_edge [iedg] == not_marked) continue;
marked_edge [iedg] += gamma_first_dis_inod;
}
// 2.4) append orphan node to regular one in gamma_node_list and set marked_node
for (index_set::const_iterator
iter = orphan_node[0].begin(),
last = orphan_node[0].end(); iter != last; ++iter) {
size_type dis_inod = *iter;
marked_node.dis_entry(dis_inod) = gamma_first_dis_inod + gamma_node_list.size();
gamma_node_list.push_back (lambda.dis_node(dis_inod));
}
marked_node.dis_entry_assembly();
marked_edge.dis_entry_assembly();
// ------------------------------------------------------------
// 3) convert lists to fixed size distributed arrays
// ------------------------------------------------------------
disarray<point_basic<T>, M> gamma_node;
std::array<disarray<element_type,M>, reference_element::max_variant> gamma_side;
gamma_list2disarray (gamma_node_list, gamma_side_list, comm, d, gamma_node, gamma_side, sid_ie2bnd_ie);
// ------------------------------------------------------------
// 4) replace inod to dis_inod in element lists
// ------------------------------------------------------------
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
for (size_type ige = 0, nge = gamma_side[variant].size(); ige < nge; ige++) {
element_type& S = gamma_side[variant][ige];
for (size_type loc_inod = 0, loc_nnod = S.n_node(); loc_inod < loc_nnod; loc_inod++) {
if (S[loc_inod] == undef) continue; // external node, will be solved later
S[loc_inod] += gamma_first_dis_inod;
}
}
}
// ----------------------------------------------------------------
// 5) solve intersection that are located on external edges & nodes
// ----------------------------------------------------------------
marked_node.set_dis_indexes (ext_marked_node_idx);
for (list<to_solve>::const_iterator iter = node_to_solve.begin(),
last = node_to_solve.end(); iter != last; iter++) {
const to_solve& x = *iter;
element_type& S = gamma_side[x.variant][x.S_ige];
check_macro (S[x.k] == undef, "external index already solved");
size_type dis_inod = x.dis_i;
size_type iproc = node_ownership.find_owner(dis_inod);
size_type gamma_dis_inod = marked_node.dis_at (dis_inod);
S[x.k] = gamma_dis_inod;
}
marked_edge.set_dis_indexes (ext_marked_edge_idx);
for (list<to_solve>::const_iterator iter = edge_to_solve.begin(),
last = edge_to_solve.end(); iter != last; iter++) {
const to_solve& x = *iter;
element_type& S = gamma_side[x.variant][x.S_ige];
check_macro (S[x.k] == undef, "external index already solved");
size_type dis_iedg = x.dis_i;
size_type iproc = edge_ownership.find_owner(dis_iedg);
size_type gamma_dis_inod = marked_edge.dis_at (dis_iedg);
S[x.k] = gamma_dis_inod;
}
// ------------------------------------------------------------
// 6) convert intersection to geo
// ------------------------------------------------------------
geo_basic<T,M> gamma (lambda, gamma_node, gamma_side);
return gamma;
}
template <class T, class M>
geo_basic<T,M>
level_set (
const field_basic<T,M>& fh,
const level_set_option& opt)
{
typedef geo_element::size_type size_type;
std::vector<size_type> bnd_dom_ie_list;
disarray<size_type,M> sid_ie2bnd_ie;
return level_set_internal (fh, opt, bnd_dom_ie_list, sid_ie2bnd_ie);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M) \
template geo_basic<T,M> level_set_internal ( \
const field_basic<T,M>&, \
const level_set_option&, \
std::vector<size_t>&, \
disarray<size_t,M>&); \
template geo_basic<T,M> level_set ( \
const field_basic<T,M>&, \
const level_set_option&);
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
} // namespace
size_type
field::size_type size_type
Definition: branch.cc:425
Float
see the Float page for the full documentation
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
rheolef::disarray::iterator
rep::base::iterator iterator
Definition: disarray.h:464
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::distributor::find_owner
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
Definition: distributor.cc:106
rheolef::distributor::size
size_type size(size_type iproc) const
Definition: distributor.h:163
rheolef::distributor::first_index
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition: distributor.h:151
rheolef::distributor::is_owned
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
Definition: distributor.h:213
rheolef::distributor::decide
static const size_type decide
Definition: distributor.h:76
rheolef::distributor::comm
const communicator_type & comm() const
Definition: distributor.h:145
rheolef::field_basic< T, M >
rheolef::field_basic::get_space
const space_type & get_space() const
Definition: field.h:300
rheolef::field_basic::dis_dof
const T & dis_dof(size_type dis_idof) const
Definition: field.cc:376
rheolef::field_basic::get_approx
std::string get_approx() const
Definition: field.h:303
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition: field.h:301
rheolef::field_basic::dis_dof_update
void dis_dof_update() const
Definition: field.cc:410
rheolef::geo_basic< T, M >
rheolef::geo_element_auto
Definition: geo_element.h:376
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_element::edge
size_type edge(size_type i) const
Definition: geo_element.h:209
rheolef::geo_element::face
size_type face(size_type i) const
Definition: geo_element.h:210
rheolef::geo_element::size_type
reference_element::size_type size_type
Definition: geo_element.h:125
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::geo_element::name
char name() const
Definition: geo_element.h:169
rheolef::heap_allocator
Definition: heap_allocator.h:69
rheolef::index_set
Definition: index_set_header.icc:58
rheolef::index_set::insert
void insert(size_type dis_i)
Definition: index_set_header.icc:158
rheolef::point_basic
Definition: point.h:87
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element::max_variant
static const variant_type max_variant
Definition: reference_element.h:82
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::space_basic< T, M >
rheolef::std
Definition: vec_expr_v2.h:402
rheolef::element_type
geo_element_auto< heap_allocator< geo_element::size_type > > element_type
Definition: level_set.cc:346
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::distributed
distributed
Definition: asr.cc:228
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
rheolef::level_set_epsilon
static Float level_set_epsilon
Definition: level_set.cc:81
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
T
Expr1::float_type T
Definition: field_expr.h:261
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.f
int f
Definition: mkgeo_ball.sh:221
mkgeo_ball.variant
int variant
Definition: mkgeo_ball.sh:149
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
mkgeo_ball.a
int a
Definition: mkgeo_ball.sh:151
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef::io::out
@ out
Definition: rheostream.h:167
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
rheolef::space_numbering::dis_idof
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
Definition: space_numbering.cc:187
rheolef::space_numbering::dis_inod
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
Definition: space_numbering.cc:177
rheolef::space_numbering::nnod
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:54
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::intersection_is_quadrilateral_T
bool intersection_is_quadrilateral_T(const std::vector< T > &f, quadruplet &q)
Definition: level_set.cc:221
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
rheolef::vect
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:265
rheolef::gamma_list2disarray
void gamma_list2disarray(const std::list< point_basic< T > > &gamma_node_list, std::array< std::list< std::pair< element_type, size_t > >, reference_element::max_variant > gamma_side_list, const communicator &comm, size_t d, disarray< point_basic< T >, M > &gamma_node, std::array< disarray< element_type, M >, reference_element::max_variant > &gamma_side, disarray< size_t, M > &sid_ie2bnd_ie)
Definition: level_set.cc:350
rheolef::operator<<
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition: catchmark.h:99
rheolef::level_set_internal
geo_basic< T, M > level_set_internal(const field_basic< T, M > &, const level_set_option &, std::vector< size_t > &, disarray< size_t, M > &)
Definition: level_set.cc:415
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::level_set
geo_basic< T, M > level_set(const field_basic< T, M > &fh, const level_set_option &opt)
Definition: level_set.cc:924
rk::gamma
Float gamma[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:70
a
Definition: diffusion_isotropic.h:25
f
Definition: cavity_dg.h:29
lambda
Definition: yield_slip_circle.h:34
p
Definition: sphere.icc:25
rheolef::const_iterator
Definition: field_expr_recursive.h:552
rheolef::level_set_option
Definition: level_set.h:66
rheolef::level_set_option::epsilon
Float epsilon
Definition: level_set.h:68
rheolef::level_set_option::split_to_triangle
bool split_to_triangle
Definition: level_set.h:67
epsilon
Float epsilon
Definition: transmission_error.cc:25
M
Expr1::memory_type M
Definition: vec_expr_v2.h:416