You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
222 lines
7.5 KiB
222 lines
7.5 KiB
3 years ago
|
#define CGAL_USE_BASIC_VIEWER
|
||
|
|
|
||
|
|
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||
|
|
|
||
|
|
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
|
||
|
|
#include <CGAL/Triangulation_face_base_with_info_2.h>
|
||
|
|
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
|
||
|
|
|
||
|
|
#include <CGAL/Delaunay_mesher_2.h>
|
||
|
|
#include <CGAL/Delaunay_mesh_face_base_2.h>
|
||
|
|
#include <CGAL/Delaunay_mesh_vertex_base_2.h>
|
||
|
|
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
|
||
|
|
#include <CGAL/lloyd_optimize_mesh_2.h>
|
||
|
|
#include <CGAL/Polygon_2.h>
|
||
|
|
|
||
|
|
#include <iostream>
|
||
|
|
#include <cstdint>
|
||
|
|
#include <fstream>
|
||
|
|
#include <iostream>
|
||
|
|
|
||
|
|
#include "mex.hpp"
|
||
|
|
#include "mexAdapter.hpp"
|
||
|
|
#include "MatlabDataArray.hpp"
|
||
|
|
#include <omp.h>
|
||
|
|
|
||
|
|
struct FaceInfo2
|
||
|
|
{
|
||
|
|
FaceInfo2() {}
|
||
|
|
int nesting_level;
|
||
|
|
bool in_domain() {
|
||
|
|
return nesting_level % 2 == 1;
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||
|
|
|
||
|
|
typedef K::Point_2 Point;
|
||
|
|
//typedef CGAL::Triangulation_vertex_base_2<K> Vb1;
|
||
|
|
typedef CGAL::Triangulation_vertex_base_with_info_2<unsigned int, K> Vb1;
|
||
|
|
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K> Fbb;
|
||
|
|
typedef CGAL::Constrained_triangulation_face_base_2<K, Fbb> Fb1;
|
||
|
|
typedef CGAL::Triangulation_data_structure_2<Vb1, Fb1> TDS1;
|
||
|
|
typedef CGAL::Exact_predicates_tag Itag;
|
||
|
|
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS1, Itag> CDT1;
|
||
|
|
|
||
|
|
typedef CGAL::Delaunay_mesh_vertex_base_2<K> Vb;
|
||
|
|
typedef CGAL::Delaunay_mesh_face_base_2<K> Fb;
|
||
|
|
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> TDS;
|
||
|
|
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS> CDT;
|
||
|
|
|
||
|
|
typedef CDT1::Point Point;
|
||
|
|
typedef CGAL::Polygon_2<K> Polygon2;
|
||
|
|
typedef CDT1::Face_handle Face_handle;
|
||
|
|
|
||
|
|
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria;
|
||
|
|
typedef CGAL::Delaunay_mesher_2<CDT, Criteria> Mesher;
|
||
|
|
|
||
|
|
void mark_domains(CDT1& ct, Face_handle start, int index, std::list<CDT1::Edge>& border)
|
||
|
|
{
|
||
|
|
if (start->info().nesting_level != -1) {
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
std::list<Face_handle> queue;
|
||
|
|
queue.push_back(start);
|
||
|
|
while (!queue.empty()) {
|
||
|
|
Face_handle fh = queue.front();
|
||
|
|
queue.pop_front();
|
||
|
|
if (fh->info().nesting_level == -1) {
|
||
|
|
fh->info().nesting_level = index;
|
||
|
|
for (int i = 0; i < 3; i++) {
|
||
|
|
CDT1::Edge e(fh, i);
|
||
|
|
Face_handle n = fh->neighbor(i);
|
||
|
|
if (n->info().nesting_level == -1) {
|
||
|
|
if (ct.is_constrained(e)) border.push_back(e);
|
||
|
|
else queue.push_back(n);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
//explore set of facets connected with non constrained edges,
|
||
|
|
//and attribute to each such set a nesting level.
|
||
|
|
//We start from facets incident to the infinite vertex, with a nesting
|
||
|
|
//level of 0. Then we recursively consider the non-explored facets incident
|
||
|
|
//to constrained edges bounding the former set and increase the nesting level by 1.
|
||
|
|
//Facets in the domain are those with an odd nesting level.
|
||
|
|
void mark_domains(CDT1& cdt)
|
||
|
|
{
|
||
|
|
for (CDT1::Face_handle f : cdt.all_face_handles()) {
|
||
|
|
f->info().nesting_level = -1;
|
||
|
|
}
|
||
|
|
std::list<CDT1::Edge> border;
|
||
|
|
mark_domains(cdt, cdt.infinite_face(), 0, border);
|
||
|
|
while (!border.empty()) {
|
||
|
|
CDT1::Edge e = border.front();
|
||
|
|
border.pop_front();
|
||
|
|
Face_handle n = e.first->neighbor(e.second);
|
||
|
|
if (n->info().nesting_level == -1) {
|
||
|
|
mark_domains(cdt, n, e.first->info().nesting_level + 1, border);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
class MexFunction : public matlab::mex::Function {
|
||
|
|
public:
|
||
|
|
void operator()(matlab::mex::ArgumentList outputs, matlab::mex::ArgumentList inputs)
|
||
|
|
{
|
||
|
|
checkArguments(outputs, inputs);
|
||
|
|
|
||
|
|
matlab::data::TypedArray<double> pnts = std::move(inputs[0]);
|
||
|
|
double x = inputs[1][0];
|
||
|
|
|
||
|
|
matlab::data::CellArray nodeMat = mexDT(pnts, x);
|
||
|
|
outputs[0] = nodeMat;
|
||
|
|
}
|
||
|
|
|
||
|
|
void checkArguments(matlab::mex::ArgumentList outputs, matlab::mex::ArgumentList inputs) {
|
||
|
|
std::shared_ptr<matlab::engine::MATLABEngine> matlabPtr = getEngine();
|
||
|
|
matlab::data::ArrayFactory factory;
|
||
|
|
|
||
|
|
if (inputs.size() != 2) {
|
||
|
|
matlabPtr->feval(u"error",
|
||
|
|
0, std::vector<matlab::data::Array>({ factory.createScalar("Four inputs required") }));
|
||
|
|
}
|
||
|
|
|
||
|
|
if (inputs[0].getDimensions()[1] != 2) {
|
||
|
|
matlabPtr->feval(u"error",
|
||
|
|
0, std::vector<matlab::data::Array>({ factory.createScalar("Points input 1 must be NX2 dimension") }));
|
||
|
|
}
|
||
|
|
if (inputs[0].getType() != matlab::data::ArrayType::DOUBLE ||
|
||
|
|
inputs[0].getType() == matlab::data::ArrayType::COMPLEX_DOUBLE) {
|
||
|
|
matlabPtr->feval(u"error",
|
||
|
|
0, std::vector<matlab::data::Array>({ factory.createScalar("Points Input 1 must be a noncomplex scalar double") }));
|
||
|
|
}
|
||
|
|
|
||
|
|
if (inputs[1].getDimensions()[0] != 1 || inputs[1].getDimensions()[1] != 1) {
|
||
|
|
matlabPtr->feval(u"error", 0,
|
||
|
|
std::vector<matlab::data::Array>({ factory.createScalar("input 2 must be 1X1 dimension") }));
|
||
|
|
}
|
||
|
|
if (inputs[1].getType() != matlab::data::ArrayType::DOUBLE ||
|
||
|
|
inputs[1].getType() == matlab::data::ArrayType::COMPLEX_DOUBLE) {
|
||
|
|
matlabPtr->feval(u"error",
|
||
|
|
0, std::vector<matlab::data::Array>({ factory.createScalar("input 2 must be a noncomplex scalar double") }));
|
||
|
|
}
|
||
|
|
if (outputs.size() > 1) {
|
||
|
|
matlabPtr->feval(u"error", 0,
|
||
|
|
std::vector<matlab::data::Array>({ factory.createScalar("Only one output is returned") }));
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
matlab::data::CellArray mexDT(matlab::data::TypedArray<double> pnts, double x)
|
||
|
|
{
|
||
|
|
// step.0 generate polygon
|
||
|
|
Polygon2 pgon;
|
||
|
|
for (int i = 0; i < pnts.getDimensions()[0]; ++i)
|
||
|
|
pgon.push_back(Point(pnts[i][0], pnts[i][1]));
|
||
|
|
|
||
|
|
// step.1 lloyd opt
|
||
|
|
CDT cdt;
|
||
|
|
cdt.insert_constraint(pgon.vertices_begin(), pgon.vertices_end(), true);
|
||
|
|
|
||
|
|
Mesher mesher(cdt);
|
||
|
|
mesher.set_criteria(Criteria(0.125, x));
|
||
|
|
mesher.refine_mesh();
|
||
|
|
|
||
|
|
CGAL::lloyd_optimize_mesh_2(cdt, CGAL::parameters::max_iteration_number = 20);
|
||
|
|
|
||
|
|
// step.2 transform to CDT1
|
||
|
|
std::vector<std::pair<Point, unsigned>> points;
|
||
|
|
int index = 0;
|
||
|
|
for (CDT::Vertex_handle vh : cdt.finite_vertex_handles())
|
||
|
|
points.push_back(std::make_pair(vh->point(), index++));
|
||
|
|
|
||
|
|
CDT1 cdt1;
|
||
|
|
cdt1.insert(points.begin(), points.end());
|
||
|
|
cdt1.insert_constraint(pgon.vertices_begin(), pgon.vertices_end(), true);
|
||
|
|
mark_domains(cdt1);
|
||
|
|
|
||
|
|
// step.3 write
|
||
|
|
Eigen::MatrixXd nodes = Eigen::MatrixXd(cdt1.number_of_vertices(), 2);
|
||
|
|
std::vector<std::vector<int>> faces;
|
||
|
|
|
||
|
|
for (CDT1::Finite_faces_iterator f = cdt1.finite_faces_begin(); f != cdt1.finite_faces_end(); f++)
|
||
|
|
{
|
||
|
|
// find the faces in polygon
|
||
|
|
if (f->info().in_domain())
|
||
|
|
{
|
||
|
|
std::vector<int> face_i{0,0,0};
|
||
|
|
for (int i = 0; i < 3; ++i)
|
||
|
|
{
|
||
|
|
Point p = f->vertex(i)->point();
|
||
|
|
int idx = f->vertex(i)->info();
|
||
|
|
nodes(idx, 0) = p.x();
|
||
|
|
nodes(idx, 1) = p.y();
|
||
|
|
face_i[i] = idx;
|
||
|
|
}
|
||
|
|
faces.push_back(face_i);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
matlab::data::ArrayFactory factory;
|
||
|
|
matlab::data::TypedArray<double> nodesMat = factory.createArray<double>({ cdt1.number_of_vertices(), 2 });
|
||
|
|
for (int i = 0; i < nodes.rows(); i++)
|
||
|
|
{
|
||
|
|
nodesMat[i][0] = nodes(i, 0);
|
||
|
|
nodesMat[i][1] = nodes(i, 1);
|
||
|
|
}
|
||
|
|
|
||
|
|
matlab::data::TypedArray<double> facesMat = factory.createArray<double>({ faces.size(), 3 });
|
||
|
|
for (int i = 0; i < faces.size(); ++i)
|
||
|
|
{
|
||
|
|
facesMat[i][0] = faces[i][0] + 1;
|
||
|
|
facesMat[i][1] = faces[i][1] + 1;
|
||
|
|
facesMat[i][2] = faces[i][2] + 1;
|
||
|
|
}
|
||
|
|
|
||
|
|
matlab::data::CellArray trimesh = factory.createCellArray({ 2, 1 });
|
||
|
|
trimesh[0] = nodesMat;
|
||
|
|
trimesh[1] = facesMat;
|
||
|
|
return trimesh;
|
||
|
|
}
|
||
|
|
};
|