mkgeo_ball_gmsh_fix.cc
Go to the documentation of this file.
// fix the blending function for curved elements in gmsh:
// face boundary interior points are incorrects
// TODO:
// - fix_full : filled circle and sphere meshes
// - fix_s : when meshes by quadrangles
// - fix_s : when input geometry is an ellipse
// - support mpirun : node are distributed
//
#include "rheolef.h"
#include "rheolef/geo_element_curved_ball.h"
#include <cstring>
using namespace rheolef;
point project_on_unit_sphere (const point& x) { return x/norm(x); }
geo_basic<Float,sequential> fix_s (const geo_basic<Float,sequential>& gamma)
{
typedef geo_basic<Float,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
size_t order = gamma.order();
disarray<point,sequential> node = gamma.get_nodes();
// 1) put all nodes exactly on the unit sphere
for (size_type iv = 0, nv = gamma.n_vertex(); iv < nv; iv++) {
node[iv] = project_on_unit_sphere (node[iv]);
}
// 2) fix nodes along edges:
std::vector<size_type> dis_inod;
for (size_type iedg = 0, nedg = gamma.size(1); iedg < nedg; iedg++) {
const geo_element& E = gamma.get_geo_element (1, iedg);
gamma.dis_inod (E, dis_inod);
const point& x0 = node[dis_inod[0]];
const point& x1 = node[dis_inod[1]];
for (size_type i = 1; i < order; i++) {
size_type loc_inod = reference_element_e::ilat2loc_inod (order, ilat(i));
point xi = x0 + ((1.0*i)/order)*(x1 - x0);
node[dis_inod[loc_inod]] = project_on_unit_sphere (xi);
}
}
// 3) fix nodes inside elements
for (size_type ie = 0, ne = gamma.size(2); ie < ne; ie++) {
const geo_element& K = gamma.get_geo_element (2, ie);
gamma.dis_inod (K, dis_inod);
switch (K.variant()) {
case reference_element::t: {
const point& x0 = node[dis_inod[0]];
const point& x1 = node[dis_inod[1]];
const point& x2 = node[dis_inod[2]];
for (size_type i = 1; i < order; i++) {
for (size_type j = 1; i+j < order; j++) {
size_type loc_inod = reference_element_t::ilat2loc_inod (order, ilat(i, j));
point hat_x ((1.0*i)/order, (1.0*j)/order);
point xi = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
node[dis_inod[loc_inod]] = project_on_unit_sphere (xi);
}
}
break;
}
case reference_element::q: {
const point& x0 = node[dis_inod[0]];
const point& x1 = node[dis_inod[1]];
const point& x2 = node[dis_inod[2]];
const point& x3 = node[dis_inod[3]];
for (size_type i = 1; i < order; i++) {
for (size_type j = 1; j < order; j++) {
size_type loc_inod = reference_element_q::ilat2loc_inod (order, ilat(i, j));
point hat_x (-1+2.*i/order, -1+2.*j/order);
point xi = 0.25*((1 - hat_x[0])*(1 - hat_x[1])*x0
+ (1 + hat_x[0])*(1 - hat_x[1])*x1
+ (1 + hat_x[0])*(1 + hat_x[1])*x2
+ (1 - hat_x[0])*(1 + hat_x[1])*x3);
node[dis_inod[loc_inod]] = project_on_unit_sphere (xi);
}
}
break;
}
default: error_macro ("element variant `"<<K.name()<<"' not yet supported");
}
}
geo_basic<Float,sequential> gamma_out = gamma;
gamma_out.set_nodes (node);
return gamma_out;
}
geo_basic<Float,sequential> fix_filled (const geo_basic<Float,sequential>& omega) {
typedef geo_basic<Float,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
static const Float pi = acos(Float(-1.0));
size_t order = omega.order();
// 1) loop on the boundary and mark boundary faces, edges & vertices
const geo_basic<Float,sequential>& bdry = omega["boundary"];
const size_type d = omega.map_dimension();
const size_type sid_dim = bdry.map_dimension();
check_macro (d == sid_dim+1, "unexpected dimensions");
std::vector<bool> bdry_ver_mark (omega.n_node(), false);
std::vector<bool> bdry_edg_mark (omega.size(1), false);
std::vector<bool> bdry_sid_mark (omega.size(sid_dim), false);
for (size_type ioisid = 0, noisid = bdry.size(); ioisid < noisid; ioisid++) {
const geo_element& S = bdry.get_geo_element(sid_dim, ioisid);
bdry_sid_mark [S.dis_ie()] = true;
for (size_type loc_iv = 0, loc_nv = S.size(); loc_iv < loc_nv; loc_iv++) {
bdry_ver_mark [S[loc_iv]] = true;
}
if (d != 3) continue;
for (size_type loc_iedg = 0, loc_nedg = S.n_subgeo(1); loc_iedg < loc_nedg; loc_iedg++) {
bdry_edg_mark [S.edge(loc_iedg)] = true;
}
}
// 2) fix boundary vertices
disarray<point,sequential> node = omega.get_nodes();
for (size_type iv = 0, nv = omega.n_vertex(); iv < nv; iv++) {
if (! bdry_ver_mark [iv]) continue;
node [iv] = project_on_unit_sphere (node[iv]);
}
// 3) fix boundary edges
std::vector<size_type> dis_inod;
for (size_type iedg = 0, nedg = omega.size(1); iedg < nedg; iedg++) {
const geo_element& E = omega.get_geo_element(1, iedg);
omega.dis_inod (E, dis_inod);
const point& x0 = node[dis_inod[0]];
const point& x1 = node[dis_inod[1]];
for (size_type i = 1; i < order; i++) {
size_type loc_inod = reference_element_e::ilat2loc_inod (order, ilat(i));
point xi = x0 + ((1.0*i)/order)*(x1 - x0);
if ((d == 2 && bdry_sid_mark [iedg]) || (d == 3 && bdry_edg_mark [iedg])) {
xi = project_on_unit_sphere (xi);
}
node[dis_inod[loc_inod]] = xi;
}
}
// 4) fix boundary sides and internal sides with one boundary edge
if (d == 3) {
for (size_type ifac = 0, nfac = omega.size(2); ifac < nfac; ifac++) {
const geo_element& S = omega.get_geo_element(2, ifac);
omega.dis_inod (S, dis_inod);
bool side_has_boundary_edge = false;
size_type loc_bdry_iedg = 0;
if (! bdry_sid_mark [ifac]) {
size_type n_bdry_edg = 0;
for (size_type loc_iedg = 0, loc_nedg = S.n_subgeo (1); loc_iedg < loc_nedg; loc_iedg++) {
size_type iedg = S.edge(loc_iedg);
if (bdry_edg_mark [iedg]) {
n_bdry_edg++;
loc_bdry_iedg = loc_iedg;
side_has_boundary_edge = true;
}
}
check_macro (n_bdry_edg <= 1, "unsupported side with "<<n_bdry_edg<<" boundary edges");
}
switch (S.variant()) {
case reference_element::t: {
// side is internal and have a curved edge
const point& x0 = node[dis_inod[0]];
const point& x1 = node[dis_inod[1]];
const point& x2 = node[dis_inod[2]];
curved_ball_t<Float> F (x0, x1, x2, loc_bdry_iedg);
for (size_type i = 1; i < order; i++) {
for (size_type j = 1; i+j < order; j++) {
size_type loc_inod = reference_element_t::ilat2loc_inod (order, ilat(i, j));
point hat_x ((1.0*i)/order, (1.0*j)/order);
point xi;
if (side_has_boundary_edge) {
xi = F (hat_x);
} else {
xi = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
if (bdry_sid_mark [ifac]) {
xi = project_on_unit_sphere (xi);
}
}
node[dis_inod[loc_inod]] = xi;
}
}
break;
}
#ifdef TODO
case reference_element::q: {
break;
}
#endif // TODO
default: error_macro ("element variant `"<<S.name()<<"' not yet supported");
}
}
}
// 5) fix interior nodes of elements that have a curved boundary side
for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
const geo_element& K = omega.get_geo_element(d, ie);
omega.dis_inod (K, dis_inod);
// 5.1) has K one side exactly on the boundary ?
size_type n_bdry_sid = 0;
size_type loc_bdry_isid = 0;
for (size_type loc_isid = 0, loc_nsid = K.n_subgeo (sid_dim); loc_isid < loc_nsid; loc_isid++) {
size_type isid = (sid_dim == 1) ? K.edge(loc_isid) : K.face(loc_isid);
if (bdry_sid_mark [isid]) {
n_bdry_sid++;
loc_bdry_isid = loc_isid;
}
}
check_macro (n_bdry_sid <= 1, "unsupported element with "<<n_bdry_sid<<" boundary sides");
bool is_on_bdry = (n_bdry_sid == 1); // n_bdry_sid == 1 i.e. K has exactly one boundary curved side
// 5.2) has K one edge exactly on the boundary ?
size_type n_bdry_edg = 0;
size_type loc_bdry_iedg = 0;
if (d == 3 && !is_on_bdry) {
for (size_type loc_iedg = 0, loc_nedg = K.n_subgeo (1); loc_iedg < loc_nedg; loc_iedg++) {
size_type iedg = K.edge(loc_iedg);
if (bdry_edg_mark [iedg]) {
n_bdry_edg++;
loc_bdry_iedg = loc_iedg;
}
}
check_macro (n_bdry_edg <= 1, "unsupported element with "<<n_bdry_edg<<" boundary edges");
}
bool has_bdry_edge = (n_bdry_edg == 1);
switch (K.variant()) {
case reference_element::t: {
// has a curved boundary edge:
const point& x0 = node[dis_inod[0]];
const point& x1 = node[dis_inod[1]];
const point& x2 = node[dis_inod[2]];
curved_ball_t<Float> F (x0, x1, x2, loc_bdry_isid);
for (size_type i = 1; i < order; i++) {
for (size_type j = 1; i+j < order; j++) {
size_type loc_inod = reference_element_t::ilat2loc_inod (order, ilat(i, j));
point hat_x ((1.0*i)/order, (1.0*j)/order);
point x;
if (is_on_bdry) {
x = F (hat_x);
} else {
x = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
}
node[dis_inod[loc_inod]] = x;
}
}
break;
}
case reference_element::T: {
if (is_on_bdry || has_bdry_edge) { // have a triangular face or an edge on the boundary
const point& x0 = node[dis_inod[0]];
const point& x1 = node[dis_inod[1]];
const point& x2 = node[dis_inod[2]];
const point& x3 = node[dis_inod[3]];
curved_ball_T<Float> F (x0, x1, x2, x3);
if (is_on_bdry) {
F.set_boundary_face (loc_bdry_isid);
} else if (has_bdry_edge) {
F.set_boundary_edge (loc_bdry_iedg);
}
for (size_type i = 1; i < order; i++) {
for (size_type j = 1; i+j < order; j++) {
for (size_type k = 1; i+j+k < order; k++) {
size_type loc_inod = reference_element_T::ilat2loc_inod (order, ilat(i, j, k));
point hat_x ((1.0*i)/order, (1.0*j)/order, (1.0*k)/order);
node[dis_inod[loc_inod]] = F (hat_x);
}
}
}
} else { // fully interior tetra:
const point& x0 = node[dis_inod[0]];
const point& x1 = node[dis_inod[1]];
const point& x2 = node[dis_inod[2]];
const point& x3 = node[dis_inod[3]];
for (size_type i = 1; i < order; i++) {
for (size_type j = 1; i+j < order; j++) {
for (size_type k = 1; i+j+k < order; k++) {
size_type loc_inod = reference_element_T::ilat2loc_inod (order, ilat(i, j, k));
point hat_x ((1.0*i)/order, (1.0*j)/order, (1.0*k)/order);
point x = (1 - hat_x[0] - hat_x[1] - hat_x[2])*x0 + hat_x[0]*x1 + hat_x[1]*x2 + hat_x[2]*x3;
node[dis_inod[loc_inod]] = x;
}
}
}
}
break;
}
case reference_element::q: {
size_type coord[4];
size_type shift[4];
shift[0] = loc_bdry_isid;
shift[1] = (loc_bdry_isid+1)%4;
shift[2] = (loc_bdry_isid+2)%4;
shift[3] = (loc_bdry_isid+3)%4;
// [x0,x1] is the curved boundary edge:
const point& x0 = node[dis_inod[shift[0]]];
const point& x1 = node[dis_inod[shift[1]]];
const point& x2 = node[dis_inod[shift[2]]];
const point& x3 = node[dis_inod[shift[3]]];
curved_ball_q<Float> F (x0, x1, x2, x3);
for (size_type i = 1; i < order; i++) {
for (size_type j = 1; j < order; j++) {
size_type loc_inod = reference_element_q::ilat2loc_inod (order, ilat(i, j));
coord [0] = i;
coord [1] = j;
coord [2] = order - i;
coord [3] = order - j;
size_type i1 = coord [shift[0]%4];
size_type j1 = coord [shift[1]%4];
point hat_x (-1+2.0*i1/order, -1+2.0*j1/order);
point x;
if (is_on_bdry) {
x = F (hat_x);
} else {
x = 0.25*( (1 - hat_x[0])*(1 - hat_x[1])*x0
+ (1 + hat_x[0])*(1 - hat_x[1])*x1
+ (1 + hat_x[0])*(1 + hat_x[1])*x2
+ (1 - hat_x[0])*(1 + hat_x[1])*x3);
}
node[dis_inod[loc_inod]] = x;
}
}
break;
}
default: error_macro ("element variant `"<<K.name()<<"' not yet supported");
}
}
geo_basic<Float,sequential> omega_out = omega;
omega_out.set_nodes (node);
return omega_out;
}
geo_basic<Float,sequential> fix (const geo_basic<Float,sequential>& omega) {
if (omega.map_dimension() < omega.dimension()) {
return fix_s (omega);
} else {
return fix_filled (omega);
}
}
void usage() {
std::cerr << "mkgeo_ball_gmsh_fix: usage:" << std::endl
<< "mkgeo_ball_gmsh_fix "
<< "[-order int] < input.geo > output.geo"
;
exit (1);
}
int main(int argc, char**argv) {
environment distributed (argc, argv);
check_macro (communicator().size() == 1,
argv[0] << ": command may be used as mono-process only");
// ----------------------------
// scan the command line
// ----------------------------
size_t order = std::numeric_limits<size_t>::max();
for (int i = 1; i < argc; i++) {
if (strcmp (argv[i], "-order") == 0) {
if (i == argc-1) { std::cerr << argv[0] << " -order: option argument missing" << std::endl; usage(); }
order = atoi(argv[++i]);
}
else { usage(); }
}
// ----------------------------
// treatment
// ----------------------------
geo_basic<Float,sequential> omega; din >> omega;
if (order != std::numeric_limits<size_t>::max()) {
omega.reset_order (order);
}
omega = fix (omega);
dout << omega;
}
size_type
field::size_type size_type
Definition: branch.cc:425
Float
see the Float page for the full documentation
point
see the point page for the full documentation
rheolef::curved_ball_T
Definition: geo_element_curved_ball.h:174
rheolef::curved_ball_T::set_boundary_face
void set_boundary_face(size_t loc_ifac_curved)
Definition: geo_element_curved_ball.h:190
rheolef::curved_ball_T::set_boundary_edge
void set_boundary_edge(size_t loc_iedg_curved)
Definition: geo_element_curved_ball.h:205
rheolef::curved_ball_q
Definition: geo_element_curved_ball.h:55
rheolef::curved_ball_t
Definition: geo_element_curved_ball.h:111
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
rheolef::environment
see the environment page for the full documentation
Definition: environment.h:115
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
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
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::geo_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
rheolef::geo_element::name
char name() const
Definition: geo_element.h:169
rheolef::point_basic
Definition: point.h:87
rheolef::reference_element_T::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:582
rheolef::reference_element_e::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:350
rheolef::reference_element_q::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:481
rheolef::reference_element_t::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:413
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
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::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::distributed
distributed
Definition: asr.cc:228
rheolef::Float
double Float
see the Float page for the full documentation
Definition: Float.h:143
rheolef::din
idiststream din
see the diststream page for the full documentation
Definition: diststream.h:427
rheolef::dout
odiststream dout(cout)
see the diststream page for the full documentation
Definition: diststream.h:430
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
fix_s
geo_basic< Float, sequential > fix_s(const geo_basic< Float, sequential > &gamma)
Definition: mkgeo_ball_gmsh_fix.cc:34
usage
void usage()
Definition: mkgeo_ball_gmsh_fix.cc:342
main
int main(int argc, char **argv)
Definition: mkgeo_ball_gmsh_fix.cc:349
project_on_unit_sphere
point project_on_unit_sphere(const point &x)
Definition: mkgeo_ball_gmsh_fix.cc:33
fix_filled
geo_basic< Float, sequential > fix_filled(const geo_basic< Float, sequential > &omega)
Definition: mkgeo_ball_gmsh_fix.cc:101
fix
geo_basic< Float, sequential > fix(const geo_basic< Float, sequential > &omega)
Definition: mkgeo_ball_gmsh_fix.cc:335
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
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
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
rk::gamma
Float gamma[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:70
rheolef.h
rheolef - reference manual