Select Git revision
gmshRegion.cpp
-
Christophe Geuzaine authoredChristophe Geuzaine authored
mainCartesian.cpp 11.56 KiB
// lx ly lz rmax levels refcs
//
// plaqueEp.stp 0.2 0.2 0.2 0.3 3
// plaqueEpRotated.stp 0.3 0.3 0.3 0.3 3
// jonction_collee2.stp 6 6 6 10 3
// panneau_raidi_simple.stp 3 3 50 5 2 0
// plaque_trouee.stp 1 1 1 2 3
#include "Gmsh.h"
#include "GModel.h"
#include "Numeric.h"
#include "GmshMessage.h"
#include "cartesian.h"
static void insertActiveCells(double x, double y, double z, double rmax,
cartesianBox<double> &box)
{
int id1 = box.getCellContainingPoint(x - rmax, y - rmax, z - rmax);
int id2 = box.getCellContainingPoint(x + rmax, y + rmax, z + rmax);
int i1, j1 ,k1;
box.getCellIJK(id1, i1, j1, k1);
int i2, j2, k2;
box.getCellIJK(id2, i2, j2, k2);
for (int i = i1; i <= i2; i++)
for (int j = j1; j <= j2; j++)
for (int k = k1; k <= k2; k++)
box.insertActiveCell(box.getCellIndex(i, j, k));
}
static void computeLevelset(GModel *gm, cartesianBox<double> &box)
{
// tolerance for desambiguation
const double tol = box.getLC() * 1.e-12;
std::vector<SPoint3> nodes;
std::vector<int> indices;
for (cartesianBox<double>::valIter it = box.nodalValuesBegin();
it != box.nodalValuesEnd(); ++it){
nodes.push_back(box.getNodeCoordinates(it->first));
indices.push_back(it->first);
}
Msg::Info(" %d nodes in the grid at level %d", (int)nodes.size(), box.getLevel());
std::vector<double> dist, localdist;
std::vector<SPoint3> dummy;
for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
for (int i = 0; i < (*fit)->stl_triangles.size(); i += 3){
int i1 = (*fit)->stl_triangles[i];
int i2 = (*fit)->stl_triangles[i + 1];
int i3 = (*fit)->stl_triangles[i + 2];
GPoint p1 = (*fit)->point((*fit)->stl_vertices[i1]);
GPoint p2 = (*fit)->point((*fit)->stl_vertices[i2]);
GPoint p3 = (*fit)->point((*fit)->stl_vertices[i3]);
SPoint2 p = ((*fit)->stl_vertices[i1] + (*fit)->stl_vertices[i2] +
(*fit)->stl_vertices[i3]) * 0.33333333;
SVector3 N = (*fit)->normal(p);
SPoint3 P1(p1.x(), p1.y(), p1.z());
SPoint3 P2(p2.x(), p2.y(), p2.z());
SPoint3 P3(p3.x(), p3.y(), p3.z());
SVector3 NN(crossprod(P2 - P1, P3 - P1));
if (dot(NN, N) > 0)
signedDistancesPointsTriangle(localdist, dummy, nodes, P1, P2, P3);
else
signedDistancesPointsTriangle(localdist, dummy, nodes, P2, P1, P3);
if(dist.empty())
dist = localdist;
else{
for (unsigned int j = 0; j < localdist.size(); j++){
// FIXME: if there is an ambiguity assume we are inside (to
// avoid holes in the structure). This is definitely just a
// hack, as it could create pockets of matter outside the
// structure...