Skip to content
Snippets Groups Projects
Commit 8c69b58e authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

Revert "3D lloyds et rtree" 

This reverts commit e97813c4f20e55cf0d68aa2bebb8f9eef55c6822.
parent 6ea5286a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Mesh/periodical.cpp deleted 100755 → 0
+ 0
− 154
View file @ 6ea5286a
// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include "voro++.hh"
#include "periodical.h"
#include "GModel.h"
#include "meshGRegion.h"
#include <fstream>
#include <algorithm>
#include "MElement.h"
using namespace voro;
/*********class voroMetal3D*********/
voroMetal3D::voroMetal3D(){}
voroMetal3D::~voroMetal3D(){}
void voroMetal3D::execute(){
GRegion* gr;
GModel* model = GModel::current();
GModel::riter it;
for(it=model->firstRegion();it!=model->lastRegion();it++)
{
gr = *it;
if(gr->getNumMeshElements()>0){
execute(gr);
}
}
}
void voroMetal3D::execute(GRegion* gr){
int i;
int j;
MElement* element;
MVertex* vertex;
std::vector<SPoint3> vertices2;
std::set<MVertex*> vertices;
std::set<MVertex*>::iterator it;
for(i=0;i<gr->getNumMeshElements();i++){
element = gr->getMeshElement(i);
for(j=0;j<element->getNumVertices();j++){
vertex = element->getVertex(j);
vertices.insert(vertex);
}
}
for(it=vertices.begin();it!=vertices.end();it++){
vertices2.push_back(SPoint3((*it)->x(),(*it)->y(),(*it)->z()));
}
execute(vertices2);
}
void voroMetal3D::execute(std::vector<SPoint3>& vertices){
int i;
int j;
int start;
int end;
int index1;
int index2;
double x,y,z;
double x1,y1,z1;
double x2,y2,z2;
double delta;
double min_x,max_x;
double min_y,max_y;
double min_z,max_z;
voronoicell_neighbor* pointer;
voronoicell_neighbor cell;
std::vector<int> faces;
std::vector<double> voronoi_vertices;
std::vector<voronoicell_neighbor*> pointers;
std::vector<SPoint3> generators;
min_x = 1000000000.0;
max_x = -1000000000.0;
min_y = 1000000000.0;
max_y = -1000000000.0;
min_z = 1000000000.0;
max_z = -1000000000.0;
for(i=0;i<vertices.size();i++){
min_x = std::min(vertices[i].x(),min_x);
max_x = std::max(vertices[i].x(),max_x);
min_y = std::min(vertices[i].y(),min_y);
max_y = std::max(vertices[i].y(),max_y);
min_z = std::min(vertices[i].z(),min_z);
max_z = std::max(vertices[i].z(),max_z);
}
delta = 0.2*(max_x - min_x);
container cont(min_x-delta,max_x+delta,min_y-delta,max_y+delta,min_z-delta,max_z+delta,6,6,6,true,true,true,vertices.size());
for(i=0;i<vertices.size();i++){
cont.put(i,vertices[i].x(),vertices[i].y(),vertices[i].z());
}
c_loop_all loop(cont);
loop.start();
do{
cont.compute_cell(cell,loop);
loop.pos(x,y,z);
pointer = new voronoicell_neighbor();
*pointer = cell;
pointers.push_back(pointer);
generators.push_back(SPoint3(x,y,z));
}while(loop.inc());
std::ofstream file("cells.pos");
file << "View \"test\" {\n";
for(i=0;i<pointers.size();i++){
faces.clear();
voronoi_vertices.clear();
pointers[i]->face_vertices(faces);
pointers[i]->vertices(generators[i].x(),generators[i].y(),generators[i].z(),voronoi_vertices);
end = 0;
while(end<faces.size()){
start = end + 1;
end = start + faces[end];
for(j=start;j<end;j++){
if(j<end-1){
index1 = faces[j];
index2 = faces[j+1];
}
else{
index1 = faces[start];
index2 = faces[end-1];
}
x1 = voronoi_vertices[3*index1];
y1 = voronoi_vertices[3*index1+1];
z1 = voronoi_vertices[3*index1+2];
x2 = voronoi_vertices[3*index2];
y2 = voronoi_vertices[3*index2+1];
z2 = voronoi_vertices[3*index2+2];
print_segment(SPoint3(x1,y1,z1),SPoint3(x2,y2,z2),file);
}
}
}
file << "};\n";
for(i=0;i<pointers.size();i++) delete pointers[i];
}
void voroMetal3D::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file){
file << "SL ("
<< p1.x() << ", " << p1.y() << ", " << p1.z() << ", "
<< p2.x() << ", " << p2.y() << ", " << p2.z()
<< "){10, 20};\n";
}
\ No newline at end of file
Mesh/periodical.h deleted 100755 → 0
+ 0
− 16
View file @ 6ea5286a
// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include "GRegion.h"
class voroMetal3D{
public:
voroMetal3D();
~voroMetal3D();
void execute();
void execute(GRegion*);
void execute(std::vector<SPoint3>&);
void print_segment(SPoint3,SPoint3,std::ofstream&);
};
\ No newline at end of file
Mesh/simple3D.cpp deleted 100755 → 0
+ 0
− 662
View file @ 6ea5286a
// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include "simple3D.h"
#include "GModel.h"
#include "MElement.h"
#include "MElementOctree.h"
#include "meshGRegion.h"
#include <queue>
#include "rtree.h"
#include <fstream>
#define k1 0.7 //k1*h is the minimal distance between two nodes
#define k2 0.5 //k2*h is the minimal distance to the boundary
#define sqrt2 1.4143
/*********definitions*********/
class Metric{
private:
double m11,m21,m31,m12,m22,m32,m13,m23,m33;
public:
Metric();
~Metric();
void set_m11(double);
void set_m21(double);
void set_m31(double);
void set_m12(double);
void set_m22(double);
void set_m32(double);
void set_m13(double);
void set_m23(double);
void set_m33(double);
double get_m11();
double get_m21();
double get_m31();
double get_m12();
double get_m22();
double get_m32();
double get_m13();
double get_m23();
double get_m33();
};
class Node{
private:
double h;
Metric m;
SPoint3 point;
public:
double min[3];
double max[3];
Node();
Node(SPoint3);
~Node();
void set_size(double);
void set_metric(Metric);
void set_point(SPoint3);
double get_size();
Metric get_metric();
SPoint3 get_point();
};
class Wrapper{
private:
bool too_close;
Node* spawn;
Node* parent;
public:
Wrapper();
Wrapper(Node*,Node*);
~Wrapper();
void set_too_close(bool);
void set_spawn(Node*);
void set_parent(Node*);
bool get_too_close();
Node* get_spawn();
Node* get_parent();
};
/*********functions*********/
double max(double a,double b){
if(a>b) return a;
else return b;
}
double uniform_distance(SPoint3 p1,SPoint3 p2,Metric m){
double det;
double distance;
double a,b,c,d,e,f,g,h,i;
double x1,y1,z1;
double x2,y2,z2;
det = m.get_m11()*(m.get_m22()*m.get_m33() - m.get_m32()*m.get_m23());
det = det - m.get_m12()*(m.get_m21()*m.get_m33() - m.get_m31()*m.get_m23());
det = det + m.get_m13()*(m.get_m21()*m.get_m32() - m.get_m31()*m.get_m22());
a = (m.get_m22()*m.get_m33() - m.get_m32()*m.get_m23())/det;
b = (m.get_m13()*m.get_m32() - m.get_m33()*m.get_m12())/det;
c = (m.get_m12()*m.get_m23() - m.get_m22()*m.get_m13())/det;
d = (m.get_m23()*m.get_m31() - m.get_m33()*m.get_m21())/det;
e = (m.get_m11()*m.get_m33() - m.get_m31()*m.get_m13())/det;
f = (m.get_m13()*m.get_m21() - m.get_m23()*m.get_m11())/det;
g = (m.get_m21()*m.get_m32() - m.get_m31()*m.get_m22())/det;
h = (m.get_m12()*m.get_m31() - m.get_m32()*m.get_m11())/det;
i = (m.get_m11()*m.get_m22() - m.get_m21()*m.get_m12())/det;
x1 = a*p1.x() + b*p1.y() + c*p1.z();
y1 = d*p1.x() + e*p1.y() + f*p1.z();
z1 = g*p1.x() + h*p1.y() + i*p1.z();
x2 = a*p2.x() + b*p2.y() + c*p2.z();
y2 = d*p2.x() + e*p2.y() + f*p2.z();
z2 = g*p2.x() + h*p2.y() + i*p2.z();
distance = max(max(fabs(x2-x1),fabs(y2-y1)),fabs(z2-z1));
return distance;
}
bool rtree_callback(Node* neighbour,void* w){
double h;
double distance;
Metric m;
Node *spawn,*parent;
Wrapper* wrapper;
wrapper = static_cast<Wrapper*>(w);
spawn = wrapper->get_spawn();
parent = wrapper->get_parent();
h = spawn->get_size();
m = spawn->get_metric();
if(neighbour!=parent){
distance = uniform_distance(spawn->get_point(),neighbour->get_point(),m);
if(distance<k1*h){
wrapper->set_too_close(1);
return false;
}
}
return true;
}
/*********class Metric*********/
Metric::Metric(){
m11 = 1.0;
m21 = 0.0;
m31 = 0.0;
m12 = 0.0;
m22 = 1.0;
m32 = 0.0;
m13 = 0.0;
m23 = 0.0;
m33 = 1.0;
}
Metric::~Metric(){}
void Metric::set_m11(double new_m11){
m11 = new_m11;
}
void Metric::set_m21(double new_m21){
m21 = new_m21;
}
void Metric::set_m31(double new_m31){
m31 = new_m31;
}
void Metric::set_m12(double new_m12){
m12 = new_m12;
}
void Metric::set_m22(double new_m22){
m22 = new_m22;
}
void Metric::set_m32(double new_m32){
m32 = new_m32;
}
void Metric::set_m13(double new_m13){
m13 = new_m13;
}
void Metric::set_m23(double new_m23){
m23 = new_m23;
}
void Metric::set_m33(double new_m33){
m33 = new_m33;
}
double Metric::get_m11(){
return m11;
}
double Metric::get_m21(){
return m21;
}
double Metric::get_m31(){
return m31;
}
double Metric::get_m12(){
return m12;
}
double Metric::get_m22(){
return m22;
}
double Metric::get_m32(){
return m32;
}
double Metric::get_m13(){
return m13;
}
double Metric::get_m23(){
return m23;
}
double Metric::get_m33(){
return m33;
}
/*********class Node*********/
Node::Node(){}
Node::Node(SPoint3 new_point){
point = new_point;
}
Node::~Node(){}
void Node::set_size(double new_h){
h = new_h;
}
void Node::set_metric(Metric new_m){
m = new_m;
}
void Node::set_point(SPoint3 new_point){
point = new_point;
}
double Node::get_size(){
return h;
}
Metric Node::get_metric(){
return m;
}
SPoint3 Node::get_point(){
return point;
}
/*********class Wrapper*********/
Wrapper::Wrapper(){
too_close = 0;
}
Wrapper::Wrapper(Node* new_spawn,Node* new_parent){
too_close = 0;
spawn = new_spawn;
parent = new_parent;
}
Wrapper::~Wrapper(){}
void Wrapper::set_too_close(bool new_too_close){
too_close = new_too_close;
}
void Wrapper::set_spawn(Node* new_spawn){
spawn = new_spawn;
}
void Wrapper::set_parent(Node* new_parent){
parent = new_parent;
}
bool Wrapper::get_too_close(){
return too_close;
}
Node* Wrapper::get_spawn(){
return spawn;
}
Node* Wrapper::get_parent(){
return parent;
}
/*********class Filler*********/
Filler::Filler(){}
Filler::~Filler(){}
void Filler::treat_model(){
GRegion* gr;
GModel* model = GModel::current();
GModel::riter it;
for(it=model->firstRegion();it!=model->lastRegion();it++)
{
gr = *it;
if(gr->getNumMeshElements()>0){
treat_region(gr);
}
}
}
void Filler::treat_region(GRegion* gr){
int i,j;
int count;
bool ok;
double x,y,z;
SPoint3 point;
Node *node,*spawn,*parent,*n1,*n2,*n3,*n4,*n5,*n6;
MVertex* vertex;
MElement* element;
MElementOctree* octree;
deMeshGRegion deleter;
Wrapper wrapper;
std::queue<Node*> fifo;
std::vector<Node*> data;
std::vector<Node*> garbage;
std::vector<MVertex*> boundary_vertices;
std::set<MVertex*> old_vertices;
std::set<MVertex*>::iterator it;
RTree<Node*,double,3,double> rtree;
octree = new MElementOctree(gr->model());
for(i=0;i<gr->getNumMeshElements();i++){
element = gr->getMeshElement(i);
for(j=0;j<element->getNumVertices();j++){
vertex = element->getVertex(j);
old_vertices.insert(vertex);
}
}
for(it=old_vertices.begin();it!=old_vertices.end();it++){
if((*it)->onWhat()->dim()<3){
boundary_vertices.push_back(*it);
}
}
std::ofstream file("nodes.pos");
file << "View \"test\" {\n";
for(i=0;i<boundary_vertices.size();i++){
x = boundary_vertices[i]->x();
y = boundary_vertices[i]->y();
z = boundary_vertices[i]->z();
node = new Node(SPoint3(x,y,z));
compute_parameters(node);
rtree.Insert(node->min,node->max,node);
fifo.push(node);
print_node(node,file);
}
file << "};\n";
count = 1;
while(!fifo.empty()){
parent = fifo.front();
fifo.pop();
garbage.push_back(parent);
n1 = new Node();
n2 = new Node();
n3 = new Node();
n4 = new Node();
n5 = new Node();
n6 = new Node();
offsprings(octree,parent,n1,n2,n3,n4,n5,n6);
data.clear();
data.push_back(n1);
data.push_back(n2);
data.push_back(n3);
data.push_back(n4);
data.push_back(n5);
data.push_back(n6);
for(i=0;i<6;i++){
ok = 0;
spawn = data[i];
point = spawn->get_point();
x = point.x();
y = point.y();
z = point.z();
if(inside_domain(octree,x,y,z)){
compute_parameters(spawn);
if(far_from_boundary(octree,spawn)){
wrapper.set_too_close(0);
wrapper.set_spawn(spawn);
wrapper.set_parent(parent);
rtree.Search(spawn->min,spawn->max,rtree_callback,&wrapper);
if(!wrapper.get_too_close()){
fifo.push(spawn);
rtree.Insert(spawn->min,spawn->max,spawn);
vertex = new MVertex(x,y,z,gr,0);
new_vertices.push_back(vertex);
ok = 1;
}
}
}
if(!ok) delete spawn;
}
printf("%d\n",count);
count++;
}
deleter(gr);
std::vector<GRegion*> regions;
regions.push_back(gr);
meshGRegion mesher(regions); //?
mesher(gr); //?
MeshDelaunayVolume(regions);
delete octree;
for(i=0;i<garbage.size();i++) delete garbage[i];
for(i=0;i<new_vertices.size();i++) delete new_vertices[i];
new_vertices.clear();
}
Metric Filler::get_metric(double x,double y,double z){
double angle;
Metric m;
angle = atan2(z,x);
m = Metric();
m.set_m11(1.0);
m.set_m21(0.0);
m.set_m31(0.0);
m.set_m12(0.0);
m.set_m22(1.0);
m.set_m32(0.0);
m.set_m13(0.0);
m.set_m23(0.0);
m.set_m33(1.0);
return m;
}
double Filler::get_size(double x,double y,double z){
return 0.25;
}
bool Filler::inside_domain(MElementOctree* octree,double x,double y,double z){
MElement* element;
element = (MElement*)octree->find(x,y,z,3,true);
if(element!=NULL) return 1;
else return 0;
}
bool Filler::far_from_boundary(MElementOctree* octree,Node* node){
double x,y,z;
double h;
SPoint3 point;
MElement *e1,*e2,*e3,*e4,*e5,*e6;
point = node->get_point();
x = point.x();
y = point.y();
z = point.z();
h = node->get_size();
e1 = (MElement*)octree->find(x+k2*h,y,z,3,true);
e2 = (MElement*)octree->find(x-k2*h,y,z,3,true);
e3 = (MElement*)octree->find(x,y+k2*h,z,3,true);
e4 = (MElement*)octree->find(x,y-k2*h,z,3,true);
e5 = (MElement*)octree->find(x,y,z+k2*h,3,true);
e6 = (MElement*)octree->find(x,y,z-k2*h,3,true);
if(e1!=NULL && e2!=NULL && e3!=NULL && e4!=NULL && e5!=NULL && e6!=NULL) return 1;
else return 0;
}
void Filler::compute_parameters(Node* node){
double x,y,z;
double h;
Metric m;
SPoint3 point;
point = node->get_point();
x = point.x();
y = point.y();
z = point.z();
m = get_metric(x,y,z);
h = get_size(x,y,z);
node->set_size(h);
node->set_metric(m);
node->min[0] = x - sqrt2*k1*h;
node->min[1] = y - sqrt2*k1*h;
node->min[2] = z - sqrt2*k1*h;
node->max[0] = x + sqrt2*k1*h;
node->max[1] = y + sqrt2*k1*h;
node->max[2] = z + sqrt2*k1*h;
}
void Filler::offsprings(MElementOctree* octree,Node* node,Node* n1,Node* n2,Node* n3,Node* n4,Node* n5,Node* n6){
double x,y,z;
double x1,y1,z1;
double x2,y2,z2;
double x3,y3,z3;
double x4,y4,z4;
double x5,y5,z5;
double x6,y6,z6;
double h;
double h1,h2,h3,h4,h5,h6;
Metric m;
SPoint3 point;
point = node->get_point();
x = point.x();
y = point.y();
z = point.z();
h = node->get_size();
m = node->get_metric();
h1 = improvement(octree,point,h,SVector3(m.get_m11(),m.get_m21(),m.get_m31()));
x1 = x + h1*m.get_m11();
y1 = y + h1*m.get_m21();
z1 = z + h1*m.get_m31();
h2 = improvement(octree,point,h,SVector3(-m.get_m11(),-m.get_m21(),-m.get_m31()));
x2 = x - h2*m.get_m11();
y2 = y - h2*m.get_m21();
z2 = z - h2*m.get_m31();
h3 = improvement(octree,point,h,SVector3(m.get_m12(),m.get_m22(),m.get_m32()));
x3 = x + h3*m.get_m12();
y3 = y + h3*m.get_m22();
z3 = z + h3*m.get_m32();
h4 = improvement(octree,point,h,SVector3(-m.get_m12(),-m.get_m22(),-m.get_m32()));
x4 = x - h4*m.get_m12();
y4 = y - h4*m.get_m22();
z4 = z - h4*m.get_m32();
h5 = improvement(octree,point,h,SVector3(m.get_m13(),m.get_m23(),m.get_m33()));
x5 = x + h5*m.get_m13();
y5 = y + h5*m.get_m23();
z5 = z + h5*m.get_m33();
h6 = improvement(octree,point,h,SVector3(-m.get_m13(),-m.get_m23(),-m.get_m33()));
x6 = x - h6*m.get_m13();
y6 = y - h6*m.get_m23();
z6 = z - h6*m.get_m33();
*n1 = Node(SPoint3(x1,y1,z1));
*n2 = Node(SPoint3(x2,y2,z2));
*n3 = Node(SPoint3(x3,y3,z3));
*n4 = Node(SPoint3(x4,y4,z4));
*n5 = Node(SPoint3(x5,y5,z5));
*n6 = Node(SPoint3(x6,y6,z6));
}
double Filler::improvement(MElementOctree* octree,SPoint3 point,double h_nearer,SVector3 direction){
double x,y,z;
double average;
double h_farther;
double coeffA,coeffB;
x = point.x() + h_nearer*direction.x();
y = point.y() + h_nearer*direction.y();
z = point.z() + h_nearer*direction.z();
if(inside_domain(octree,x,y,z)){
h_farther = get_size(x,y,z);
}
else h_farther = h_nearer;
coeffA = 1.0;
coeffB = 0.2;
if(h_farther>h_nearer){
average = coeffA*h_nearer + (1.0-coeffA)*h_farther;
}
else{
average = coeffB*h_nearer + (1.0-coeffB)*h_farther;
}
return average;
}
int Filler::get_nbr_new_vertices(){
return new_vertices.size();
}
MVertex* Filler::get_new_vertex(int i){
return new_vertices[i];
}
void Filler::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file){
file << "SL ("
<< p1.x() << ", " << p1.y() << ", " << p1.z() << ", "
<< p2.x() << ", " << p2.y() << ", " << p2.z() << ")"
<< "{10, 20};\n";
}
void Filler::print_node(Node* node,std::ofstream& file){
double x,y,z;
double x1,y1,z1;
double x2,y2,z2;
double x3,y3,z3;
double x4,y4,z4;
double x5,y5,z5;
double x6,y6,z6;
double h;
Metric m;
SPoint3 point;
point = node->get_point();
x = point.x();
y = point.y();
z = point.z();
h = node->get_size();
m = node->get_metric();
x1 = x + h*m.get_m11();
y1 = y + h*m.get_m21();
z1 = z + h*m.get_m31();
x2 = x - h*m.get_m11();
y2 = y - h*m.get_m21();
z2 = z - h*m.get_m31();
x3 = x + h*m.get_m12();
y3 = y + h*m.get_m22();
z3 = z + h*m.get_m32();
x4 = x - h*m.get_m12();
y4 = y - h*m.get_m22();
z4 = z - h*m.get_m32();
x5 = x + h*m.get_m13();
y5 = y + h*m.get_m23();
z5 = z + h*m.get_m33();
x6 = x - h*m.get_m13();
y6 = y - h*m.get_m23();
z6 = z - h*m.get_m33();
print_segment(SPoint3(x,y,z),SPoint3(x1,y1,z1),file);
print_segment(SPoint3(x,y,z),SPoint3(x2,y2,z2),file);
print_segment(SPoint3(x,y,z),SPoint3(x3,y3,z3),file);
print_segment(SPoint3(x,y,z),SPoint3(x4,y4,z4),file);
print_segment(SPoint3(x,y,z),SPoint3(x5,y5,z5),file);
print_segment(SPoint3(x,y,z),SPoint3(x6,y6,z6),file);
}
/*********static declarations*********/
std::vector<MVertex*> Filler::new_vertices;
Mesh/simple3D.h deleted 100755 → 0
+ 0
− 33
View file @ 6ea5286a
// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include "SVector3.h"
#include <list>
#include "GRegion.h"
#include "MElementOctree.h"
class Node;
class Metric;
class Filler{
private:
static std::vector<MVertex*> new_vertices;
Metric get_metric(double,double,double);
double get_size(double,double,double);
bool inside_domain(MElementOctree*,double,double,double);
bool far_from_boundary(MElementOctree*,Node*);
void compute_parameters(Node*);
void offsprings(MElementOctree*,Node*,Node*,Node*,Node*,Node*,Node*,Node*);
double improvement(MElementOctree*,SPoint3,double,SVector3);
void print_segment(SPoint3,SPoint3,std::ofstream&);
void print_node(Node*,std::ofstream&);
public:
Filler();
~Filler();
void treat_model();
void treat_region(GRegion*);
static int get_nbr_new_vertices();
static MVertex* get_new_vertex(int);
};
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment