Commit cbf2696a authored by stenger's avatar stenger
Browse files

meshconv v3.23

parent 7a10d9c2
This diff is collapsed.
......@@ -1994,6 +1994,197 @@ double BackGrid::distance_to_surface_3d(const Vertex &destination, int &element_
return min_distance;
}
//--------------------< test_projection_distances_3d >
//tests whether one of the surface elements in the segment with indices 'ix', 'iy', 'iz' has a
//distance to 'destination' which is smaller than 'min_distance'. If so, 'min_distance' is set
//to the new smallest distance.
void BackGrid::test_projection_distances_3d(const Vertex &destination, int ix, int iy, int iz,
int tid, double &min_distance, Element &el, Vertex &proj){
AugmentedList<Element*> &lst = grid[ix][iy][iz];
for(AugmentedList<Element*>::iterator it = lst.begin(),
it_end = lst.end(); it != it_end; ++it){
Element &e = **it;
if(e.operation_index[tid] != current_operation_index[tid]){
e.operation_index[tid] = current_operation_index[tid];
Vertex tmp_proj;
double distance = projection_point_triangle_3d(destination, n->vertices[e.v[0]],
n->vertices[e.v[1]], n->vertices[e.v[2]], tmp_proj);
if(distance < min_distance){
min_distance = distance;
el = e;
proj = tmp_proj;
}
}
}
}
//--------------------< projection_to_surface_3d >
double BackGrid::projection_to_surface_3d(const Vertex &destination, Element &el, Vertex &proj){
int i, x, y, z;
int sx1, sx2, sy1, sy2, sz1, sz2;
int p1, p2, q1, q2;
double min_distance = 9e100;
//determine thread-number for operation-index-disambiguation
#ifdef OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
//traverse through segments of the backgrid in concentric rectangles around 'destination'
sx1 = (int)((destination[0] - minc[0]) / seg_size[0]);
if(sx1 >= segments[0]) sx1 = segments[0] - 1;
sy1 = (int)((destination[1] - minc[1]) / seg_size[1]);
if(sy1 >= segments[1]) sy1 = segments[1] - 1;
sz1 = (int)((destination[2] - minc[2]) / seg_size[2]);
if(sz1 >= segments[2]) sz1 = segments[2] - 1;
//stop if a unit-distance is set for the segment containing 'destination'
double ud;
#ifdef OPENMP
#ifdef NO_ATOMIC
#pragma omp critical (atomic_write)
#endif
#endif
{ud = grid[sx1][sy1][sz1].unit_distance;}
if(ud > distance_threshold){
#ifdef OPENMP
#ifdef NO_ATOMIC
#pragma omp critical (atomic_write)
#endif
#endif
{el.index = grid[sx1][sy1][sz1].entry_index;}
return ud;
}
//determine maximal number of concentric rectangles around 'destination' that might have
//intersections
sx2 = segments[0] - sx1 - 1;
if(sx1 > sx2) sx2 = sx1;
sy2 = segments[1] - sy1 - 1;
if(sy1 > sy2) sy2 = sy1;
if(sx2 > sy2) sy2 = sx2;
sz2 = segments[2] - sz1 - 1;
if(sz1 > sz2) sz2 = sz1;
if(sy2 > sz2) sz2 = sy2;
//now sz2 stores the maximal number of rectangles
++current_operation_index[tid];
//check the segment containing 'destination' for nearby surface mesh elements
test_projection_distances_3d(destination, sx1, sy1, sz1, tid, min_distance, el, proj);
//successively check segments around 'destination' with increasing distance to 'destination'
for(i = 1; i <= sz2; i++){
//determine boundaries of back and front
p1 = sx1 - i;
if(p1 < 0) p1 = 0;
p2 = sx1 + i;
if(p2 >= segments[0]) p2 = segments[0] - 1;
q1 = sy1 - i;
if(q1 < 0) q1 = 0;
q2 = sy1 + i;
if(q2 >= segments[1]) q2 = segments[1] - 1;
//back
z = sz1 - i;
if(z >= 0){
for(y = q1; y <= q2; y++){
for(x = p1; x <= p2; x++) if(!grid[x][y][z].empty()){
test_projection_distances_3d(destination, x, y, z, tid, min_distance, el, proj);
}
}
}
//front
z = sz1 + i;
if(z < segments[2]){
for(y = q1; y <= q2; y++){
for(x = p1; x <= p2; x++) if(!grid[x][y][z].empty()){
test_projection_distances_3d(destination, x, y, z, tid, min_distance, el, proj);
}
}
}
//determine boundaries of bottom and top
p1 = sx1 - i;
if(p1 < 0) p1 = 0;
p2 = sx1 + i;
if(p2 >= segments[0]) p2 = segments[0] - 1;
q1 = sz1 - i + 1;
if(q1 < 0) q1 = 0;
q2 = sz1 + i - 1;
if(q2 >= segments[2]) q2 = segments[2] - 1;
//bottom
y = sy1 - i;
if(y >= 0){
for(z = q1; z <= q2; z++){
for(x = p1; x <= p2; x++) if(!grid[x][y][z].empty()){
test_projection_distances_3d(destination, x, y, z, tid, min_distance, el, proj);
}
}
}
//top
y = sy1 + i;
if(y < segments[1]){
for(z = q1; z <= q2; z++){
for(x = p1; x <= p2; x++) if(!grid[x][y][z].empty()){
test_projection_distances_3d(destination, x, y, z, tid, min_distance, el, proj);
}
}
}
//determine boundaries of left and right
p1 = sy1 - i + 1;
if(p1 < 0) p1 = 0;
p2 = sy1 + i - 1;
if(p2 >= segments[1]) p2 = segments[1] - 1;
q1 = sz1 - i + 1;
if(q1 < 0) q1 = 0;
q2 = sz1 + i - 1;
if(q2 >= segments[2]) q2 = segments[2] - 1;
//left
x = sx1 - i;
if(x >= 0){
for(z = q1; z <= q2; z++){
for(y = p1; y <= p2; y++) if(!grid[x][y][z].empty()){
test_projection_distances_3d(destination, x, y, z, tid, min_distance, el, proj);
}
}
}
//right
x = sx1 + i;
if(x < segments[0]){
for(z = q1; z <= q2; z++){
for(y = p1; y <= p2; y++) if(!grid[x][y][z].empty()){
test_projection_distances_3d(destination, x, y, z, tid, min_distance, el, proj);
}
}
}
if(min_distance < i * min_seg_size) break;
}
//set unit-distance of the segment that contains destination-vertex
if(min_distance >= distance_threshold){
#ifdef OPENMP
#ifdef NO_ATOMIC
#pragma omp critical (atomic_write)
#else
#pragma omp atomic write
#endif
#endif
grid[sx1][sy1][sz1].unit_distance = min_distance;
#ifdef OPENMP
#ifdef NO_ATOMIC
#pragma omp critical (atomic_write)
#else
#pragma omp atomic write
#endif
#endif
grid[sx1][sy1][sz1].entry_index = el.index;
}
return min_distance;
}
//--------------------< vertices_in_sphere_3d >
//returns a list of vertices lying completely inside a sphere of the specified radius around the
//specified center_vertex.
......
This diff is collapsed.
......@@ -491,6 +491,56 @@ bool tetrahedal_interpolation_test(const Vertex &p, const Vertex &a, const Verte
return true;
}
//--------------------< cross_product >
Vertex cross_product(const Vertex &a, const Vertex &b){
return Vertex(a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]);
}
//--------------------< triangular_interpolation >
double triangular_interpolation(const Vertex &p, const Vertex &a, const Vertex &b, const Vertex &c,
const double va, const double vb, const double vc){
//calculate vectors from 'p' to 'a', 'b' and 'c'
Vertex fa = a - p;
Vertex fb = b - p;
Vertex fc = c - p;
//calculate the areas and barycentric coordinates
double area_t = cross_product(a - b, a - c).length(); // main triangle area
if(abs(area_t) < epsilon){
warning("triangular_interpolation: area_t == 0");
return 0.0;
}
double area_a = cross_product(fb, fc).length(); // a's triangle area
double area_b = cross_product(fc, fa).length(); // b's triangle area
double area_c = cross_product(fa, fb).length(); // c's triangle area
//interpolate
return (va * area_a + vb * area_b + vc * area_c) / area_t;
}
//--------------------< triangular_interpolation_test >
bool triangular_interpolation_test(const Vertex &p, const Vertex &a, const Vertex &b,
const Vertex &c){
//calculate vectors from 'p' to 'a', 'b' and 'c'
Vertex fa = a - p;
Vertex fb = b - p;
Vertex fc = c - p;
//calculate the areas and barycentric coordinates
double area_t = cross_product(a - b, a - c).length(); // main triangle area
if(abs(area_t) < epsilon){
warning("triangular_interpolation: area_t == 0");
return false;
}
double area_a = cross_product(fb, fc).length(); //a's triangle area
double area_b = cross_product(fc, fa).length(); //b's triangle area
double area_c = cross_product(fa, fb).length(); //c's triangle area
//interpolate
Vertex new_p = (a * area_a + b * area_b + c * area_c) / area_t;
if(edge_length_3d(p, new_p) > 1e-7){
warning("triangular_interpolation_test: dist = " + ft(edge_length_3d(p, new_p)));
return false;
}
return true;
}
//--------------------< intersection_line_line_3d >
bool intersection_line_line_3d(const Vertex &p1, const Vertex &p2,
const Vertex &p3, const Vertex &p4){
......@@ -1094,12 +1144,12 @@ double distance_point_line_3d(const Vertex &point, const Vertex &lin0, const Ver
//if the orthogonal projection of 'point' onto l is not on the given line segment then one of the
//vertices 'lin0' or 'lin1' is the nearest point to 'point'
if(t < 0){
if(t <= 0.0){
//'lin0' is nearest
fx = lin0[0];
fy = lin0[1];
fz = lin0[2];
}else if(t > 1){
}else if(t >= 1.0){
//'lin1' is nearest
fx = lin1[0];
fy = lin1[1];
......@@ -1150,8 +1200,7 @@ double distance_point_triangle_3d(const Vertex &point,
if(abs(b1) < epsilon){
v = c2 / b2;
u = c1 / a1;
}
else{
}else{
v = (c1 * b1 - c2 * a1) / (b1 * b1 - b2 * a1);
u = (c1 - v * b1) / a1;
}
......@@ -1163,8 +1212,7 @@ double distance_point_triangle_3d(const Vertex &point,
fy = point[1] - tri0[1] - u * gy - v * hy;
fz = point[2] - tri0[2] - u * gz - v * hz;
distance = sqrt(fx * fx + fy * fy + fz * fz);
}
else{
}else{
//no, not inside. The desired distance is the distance to one of the triangle-edges
distance = 9e100;
tmp = distance_point_line_3d(point, tri0, tri1);
......@@ -1177,4 +1225,134 @@ double distance_point_triangle_3d(const Vertex &point,
return distance;
}
//--------------------< projection_point_line_3d >
//calculates the orthogonal projection from a given 'point' to the line segment between 'lin0' and
//'lin1'. If the projection lies outside the line segment the closer end-point is returned. While
//the projection is returned in argument 'proj' the actual return value of the function is the
//distance between the point and the line segment.
double projection_point_line_3d(const Vertex &point, const Vertex &lin0, const Vertex &lin1,
Vertex &proj){
//the (complete) line through 'lin0' and 'lin1' is given by 'l = lin0 + t*g' with g as follows:
Vertex g = lin1 - lin0;
//the normal-vector from 'point' to l is 'v = lin0 + t*g - point' for some t which is determined
// by 'v*g = 0'
double t = (-g[0] * (lin0[0] - point[0]) - g[1] * (lin0[1] - point[1])
- g[2] * (lin0[2] - point[2])) / (g[0] * g[0] + g[1] * g[1] + g[2] * g[2]);
//if the orthogonal projection of 'point' onto l is not on the given line segment then one of the
//vertices 'lin0' or 'lin1' is the nearest point to 'point'
if(t <= 0.0) proj = lin0; //'lin0' is nearest
else if(t >= 1.0) proj = lin1; //'lin1' is nearest
else proj = lin0 + g * t; //orthogonal projection is nearest
//calculate distance
return (point - proj).length();
}
//--------------------< projection_point_triangle_3d >
//calculates the orthogonal projection from a given 'point' to the triangle given by 'tri0',
//'tri1', 'tri2'. If the projection lies outside the triangle the closest projection to one of
//triangle's sides is returned. While the projection is returned in argument 'proj' the actual
//return value of the function is the distance between the point and the triangle.
double projection_point_triangle_3d(const Vertex &point,
const Vertex &tri0, const Vertex &tri1, const Vertex &tri2, Vertex &proj){
double gx, gy, gz, hx, hy, hz, u, v;
double c1, c2, b1, b2, a1;
//triangle lies in plane 'p = tri0 + u*g + v*h', where g and h are defined as follows:
gx = tri1[0] - tri0[0];
gy = tri1[1] - tri0[1];
gz = tri1[2] - tri0[2];
hx = tri2[0] - tri0[0];
hy = tri2[1] - tri0[1];
hz = tri2[2] - tri0[2];
//a vector from 'point' to p can be described by 'l = tri0 + u*g + v*h - point' and since we
//search the distance we need such a vector that is orthogonal to p. Thus 'l*g = 0' and 'l*h = 0'.
//We get a system of two equations of the form:
//a1*u + b1*v = c1
//b1*u + b2*v = c2
//with a1, b1, b2, c1, c2 as follows:
b1 = gx * hx + gy * hy + gz * hz;
c2 = -hx * (tri0[0] - point[0]) - hy * (tri0[1] - point[1]) - hz * (tri0[2] - point[2]);
b2 = hx * hx + hy * hy + hz * hz;
c1 = -gx * (tri0[0] - point[0]) - gy * (tri0[1] - point[1]) - gz * (tri0[2] - point[2]);
a1 = gx * gx + gy * gy + gz * gz;
//solve the system:
if(abs(b1) < epsilon){
v = c2 / b2;
u = c1 / a1;
}else{
v = (c1 * b1 - c2 * a1) / (b1 * b1 - b2 * a1);
u = (c1 - v * b1) / a1;
}
//now test whether the orthogonal projection of 'point' onto p lies inside the triangle
if(u >= 0 && v >= 0 && u + v <= 1){
//yes, inside. The length of the normal-vector is the desired distance
proj[0] = tri0[0] + u * gx + v * hx;
proj[1] = tri0[1] + u * gy + v * hy;
proj[2] = tri0[2] + u * gz + v * hz;
return (point - proj).length();
}else{
//no, not inside. The desired distance is the distance to one of the triangle-edges
double distance = 9e100;
Vertex tmp_proj;
double tmp = projection_point_line_3d(point, tri0, tri1, proj);
tmp = projection_point_line_3d(point, tri1, tri2, tmp_proj);
if(tmp < distance){
distance = tmp;
proj = tmp_proj;
}
tmp = projection_point_line_3d(point, tri2, tri0, tmp_proj);
if(tmp < distance){
distance = tmp;
proj = tmp_proj;
}
return distance;
}
}
//--------------------< make_quaternion >
//derived from https://stackoverflow.com/questions/7582398/rotate-a-vector-about-another-vector
//Note: the identity unit quaternion: quat[0] = 0.0, quat[1] = 0.0, quat[2] = 0.0, quat[3] = 1.0;
//angle is in radians.
void make_quaternion(double *quat, const Vertex &axis, double angle){
Vertex a(axis);
double rl = sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
double ca = cos(angle);
double cq = sqrt((1.0 + ca) / 2.0); //cos(acos(ca) / 2.0);
double sq = sqrt((1.0 - ca) / 2.0); //sin(acos(ca) / 2.0);
a *= sq / rl; // i.e., scaling each element, and finally:
quat[0] = a[0]; quat[1] = a[1]; quat[2] = a[2]; quat[3] = cq;
}
//--------------------< quaternion_multiplication >
//derived from https://stackoverflow.com/questions/7582398/rotate-a-vector-about-another-vector
void quaternion_multiplication(double *r, const double *q, const double *p){
double w0 = q[3], w1 = p[3];
double x0 = q[0], x1 = p[0];
double y0 = q[1], y1 = p[1];
double z0 = q[2], z1 = p[2];
r[3] = w0 * w1 - x0 * x1 - y0 * y1 - z0 * z1;
r[0] = w0 * x1 + x0 * w1 + y0 * z1 - z0 * y1;
r[1] = w0 * y1 + y0 * w1 + z0 * x1 - x0 * z1;
r[2] = w0 * z1 + z0 * w1 + x0 * y1 - y0 * x1;
}
//--------------------< quaternion_rotation >
//derived from https://stackoverflow.com/questions/7582398/rotate-a-vector-about-another-vector
//3D vector rotation: v = q * v * conj(q)
Vertex quaternion_rotation(Vertex v, const double *q){
double r[4], p[4];
r[0] = v[0]; r[1] = v[1]; r[2] = v[2]; r[3] = 0.0;
quaternion_multiplication(r, q, r);
p[0] = -q[0]; p[1] = -q[1]; p[2] = -q[2]; p[3] = q[3];
quaternion_multiplication(r, r, p);
return Vertex(r[0], r[1], r[2]);
}
} //end namespace meshconv
......@@ -156,7 +156,8 @@ Mesh* lookup_mesh(string str, bool may_create_new_mesh/*=false*/, bool no_error/
}
//--------------------< lookup_vector >
DataVector* lookup_vector(string str, bool may_create_new_vector, bool no_error){
DataVector* lookup_vector(string str, bool may_create_new_vector/*=false*/,
bool no_error/*=false*/){
map<string, DataVector >::iterator v_it;
map<string, VectorPartitioning>::iterator p_it;
......@@ -635,7 +636,7 @@ void vector_info(string name, int task_description_length){
unsigned num_partitions = vp->size();
cout << num_partitions << " partitions: " << endl;
plen = ft(num_partitions - 1).size();
for(i=0; i<num_partitions; ++i){
for(i = 0; i < num_partitions; ++i){
if((*vp)[i].second != NULL){
cout << string(task_description_length, ' ') << "partition " << ft(i, plen)
<< ": entries: " << (*vp)[i].second->size() << " (memory)" << endl;
......@@ -655,6 +656,25 @@ void vector_info(string name, int task_description_length){
}
}
//--------------------< vector_stats >
void vector_stats(string name, int task_description_length){
DataVector *v_ptr = lookup_vector(name, false, true);
if(v_ptr == NULL){
error("Vector \"" + name + "\" not declared.");
}else{
DataVector &v = *v_ptr;
double minv = infinity, maxv = minus_infinity, avg = 0.0;
for(size_t i = 0, end = v.size(); i < end; ++i){
double cv = v[i];
avg += cv;
if(cv < minv) minv = cv;
if(cv > maxv) maxv = cv;
}
avg /= v.size();
cout << "range: [" << minv << ", " << maxv << "], avg: " << avg << endl;
}
}
//--------------------< debug_output_meshes >
void debug_output_meshes(){
cout << endl << "MESHES:" << endl;
......@@ -1050,6 +1070,12 @@ void delete_all(){
map<string, MeshPartitioning>::iterator pm_it;
map<string, VectorPartitioning>::iterator vm_it;
//value-loader-files
for(pair<const string, ifstream*> &vl : value_loader_files){
vl.second->close();
delete vl.second;
}
//RefinementManagers
for(rfm_it = refinement_managers.begin(); rfm_it != refinement_managers.end(); ++rfm_it){
rfm_it->second.delete_refinement_data();
......
......@@ -421,6 +421,22 @@ void add_element(Mesh &m, const int32_t *v, const int32_t *b){//$ACL_exclude
m.elements.push_back(e);
}
//--------------------< add_element >
//adds an element to the mesh; works in worlds of arbitrary dimension; the element must be either
//a line segment, a triangle or a tetrahedron.
void add_element(Mesh &m, int v0, int v1, int v2/*=0*/, int v3/*=0*/){//$ACL_exclude
Element e;
e.rd = NULL;
e.index = m.numelements;
++m.numelements;
e.v[0] = v0;
e.v[1] = v1;
e.v[2] = v2;
e.v[3] = v3;
for(int i = 0; i < 4; ++i) e.b[i] = 0; //write_mcb_file needs it to be initialized
m.elements.push_back(e);
}
//--------------------< insert_element >
//inserts an element into the mesh right before iterator it; works in worlds of arbitrary dimension;
//the element must be either a line segment, a triangle or a tetrahedron.
......
This diff is collapsed.
......@@ -44,7 +44,7 @@
namespace meshconv{
using namespace std;
#define MESHCONV_VERSION "v3.22"
#define MESHCONV_VERSION "v3.23"
#define DEBUG_COUNT_AMBIGUOUS_RAYS 0
#define DEBUG_REFINE_BACKGRID 0
......@@ -114,9 +114,17 @@ const uint32_t compression_none = 0;
const uint32_t compression_zlib = 1;
const uint32_t compression_bzip2 = 2;
//aniso-module-transition-modes:
const int amtm_linear = 1;
const int amtm_tangens = 2;
const int amtm_accelerated = 3;
const int amtm_decelerated = 4;
const int amtm_random = 5;
//various constants
const size_t npos = string::npos;
const double pi = 3.141592653589793116;
const double pi2 = pi * 2.0;
const double invalid_value = FLT_MAX;
const size_t vector_overhead = 1000000;
const int degenerate = -666;
......@@ -138,6 +146,7 @@ extern int dijkstra_iterations;
extern string undefined_variable_warning;
extern string temp_file_folder;
extern string temp_file_prefix;
extern map<string, ifstream*> value_loader_files;
extern map<string, ElementType> element_types;
extern map<string, Mesh> meshes;
......@@ -343,7 +352,7 @@ struct DataVector{
values(size, value){/*$ACL_exclude*/}
const double& operator[](const int index) const {/*$ACL_exclude*/return values[index];}
double& operator[](const int index){/*$ACL_exclude*/return values[index];}
unsigned size(){/*$ACL_exclude*/return values.size();}
unsigned size() const{/*$ACL_exclude*/return values.size();}
void resize(unsigned size, double value=0.0){/*$ACL_exclude*/values.resize(size, value);}
void clear(){/*$ACL_exclude*/values.clear();}
void push_back(double value){/*$ACL_exclude*/
......@@ -469,6 +478,7 @@ template<typename T> struct VertexMap{
c_mid((c_begin + c_end) / 2.0), c_end(c_end), p0(NULL), p1(NULL),
p2(NULL), p3(NULL), p4(NULL), p5(NULL), p6(NULL), p7(NULL){/*$ACL_exclude*/}
VertexMap(const Mesh &m);
VertexMap(const VertexMap<T> &vm);
~VertexMap();
void clear();
void split();
......@@ -548,6 +558,7 @@ void unload_mesh(string name);
void unload_vector(string name);
void mesh_info(string name, int task_description_length=0);
void vector_info(string name, int task_description_length=0);
void vector_stats(string name, int task_description_length=0);
void delete_all();
void debug_output_meshes();
void debug_output_refinement_managers();
......@@ -600,6 +611,11 @@ double tetrahedal_interpolation(const Vertex &p, const Vertex &a, const Vertex &
const Vertex &d, const double va, const double vb, const double vc, const double vd);
bool tetrahedal_interpolation_test(const Vertex &p, const Vertex &a, const Vertex &b,
const Vertex &c, const Vertex &d);
Vertex cross_product(const Vertex &a, const Vertex &b);
double triangular_interpolation(const Vertex &p, const Vertex &a, const Vertex &b, const Vertex &c,
const double va, const double vb, const double vc);
bool triangular_interpolation_test(const Vertex &p, const Vertex &a, const Vertex &b,
const Vertex &c);
bool intersection_line_line_3d(const Vertex &p1, const Vertex &p2, const Vertex &p3,
const Vertex &p4);
bool intersection_line_line_with_intersection_3d(const Vertex &p1, const Vertex &p2,
......@@ -618,6 +634,13 @@ bool intersection_triangle_tetrahedron_3d(const Vertex &tri0, const Vertex &tri1
double distance_point_line_3d(const Vertex &point, const Vertex &lin0, const Vertex &lin1);
double distance_point_triangle_3d(const Vertex &point, const Vertex &tri0, const Vertex &tri1,
const Vertex &tri2);
double projection_point_line_3d(const Vertex &point, const Vertex &lin0, const Vertex &lin1,
Vertex &proj);
double projection_point_triangle_3d(const Vertex &point,
const Vertex &tri0, const Vertex &tri1, const Vertex &tri2, Vertex &proj);
void make_quaternion(double *quat, const Vertex &axis, double angle);
void quaternion_multiplication(double *r, const double *q, const double *p);
Vertex quaternion_rotation(Vertex v, const double *q);
//--------------------< mesh maintenance prototypes >