Skip to content

I
INMOST
Commit fb24de0b authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files
Options

Many fixes 

OldDrawGrid example: mesh slicing algorithm can now correctly resolve
multiple slices on single cell.

USE_AUTODIFF do not require USE_SOLVER and USE_MESH:
If USE_MESH is not activated, class Automatizator and algebra for
abstract variables is not compiled.
If USE_SOLVER is not activated, class Residual is not compiled.

Fixed reverse_iterator::operator- in ElementArray.

class incident_matrix that computes closed loops out of set of edges or
faces can now work on coordinates other then physical coordinates.

Fixed Face::SplitFace() was not connecting newly created faces to cells.

Cell::Volume() now uses algorithm working correctly on non-convex cells.

Face::CheckNormalOrientation() will now use volume calculation
algorithm for non-convex cells to correctly orient faces.

Face::FixNormalOrientation() will now use Face::SwapCells() instead of
Face::ReorderEdges() on internal faces. New version is thread-safe with
respect to Face::CheckNormalOrientation().

Improved algorithm for loading eclipse meshes.
parent f1025e54
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 3061 additions and 2401 deletions
Examples/OldDrawGrid/main.cpp
View file @ fb24de0b
This diff is collapsed.
Source/Autodiff/autodiff.cpp
View file @ fb24de0b
......@@ -13,6 +13,8 @@
namespace INMOST
{
#if defined(USE_MESH)
Automatizator * Automatizator::CurrentAutomatizator = NULL;
bool print_ad_ctor = false;
bool GetAutodiffPrint() {return print_ad_ctor;}
......@@ -44,6 +46,14 @@ namespace INMOST
merger.RetriveRow(r);
merger.Clear();
}
#else //USE_MESH
bool CheckCurrentAutomatizator() {return false;}
void FromBasicExpression(Sparse::Row & entries, const basic_expression & expr) {}
void AddBasicExpression(Sparse::Row & entries, INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr) {}
void FromGetJacobian(const basic_expression & expr, INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) {}
#endif //USE_MESH
#if defined(USE_MESH)
Automatizator::Automatizator(Mesh * m) :first_num(0), last_num(0), m(m) {}
Automatizator::~Automatizator()
{
......@@ -243,6 +253,7 @@ namespace INMOST
merger.Resize(0,max_unknowns,false);
#endif
}
#endif //USE_MESH
};
#endif //USE_AUTODIFF
Source/Headers/container.hpp
View file @ fb24de0b
......@@ -1757,7 +1757,7 @@ namespace INMOST
}
void push_back(const element & e)
{
{
if( pend == preserved ) reserve(capacity()*2);
new (pend++) element(e);
}
......
Source/Headers/inmost_autodiff.h
View file @ fb24de0b
......@@ -9,10 +9,10 @@
//#define NEW_VERSION
#if defined(USE_AUTODIFF) && (!defined(USE_MESH))
#warning "USE_AUTODIFF require USE_MESH"
#undef USE_AUTODIFF
#endif
//#if defined(USE_AUTODIFF) && (!defined(USE_MESH))
//#warning "USE_AUTODIFF require USE_MESH"
//#undef USE_AUTODIFF
//#endif
//#define DPRNT
......@@ -27,7 +27,7 @@ namespace INMOST
{
class Automatizator; //forward declaration
#if defined(USE_SOLVER)
class Residual
{
Sparse::Matrix jacobian;
......@@ -125,7 +125,9 @@ namespace INMOST
void UnLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.UnLock(pos);}
void TestLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.TestLock(pos);}
};
#endif
#if defined(USE_MESH)
class Automatizator
{
private:
......@@ -184,6 +186,7 @@ namespace INMOST
/// Retrive the automatizator.
static Automatizator * GetCurrent() {return CurrentAutomatizator;}
};
#endif
} //namespace INMOST
#endif //USE_AUTODIFF
......
Source/Headers/inmost_expression.h
View file @ fb24de0b
......@@ -5,11 +5,6 @@
#include <sstream> //for debug
#include <new>
#if defined(USE_AUTODIFF) && !defined(USE_SOLVER)
#warning "USE_AUTODIFF require USE_SOLVER"
#undef USE_AUTODIFF
#endif
//TODO:
// 1. Incorporate tables
......@@ -21,6 +16,10 @@
// 7. choice of directional derivatives at discontinuities for abs, pow, max, min (see ADOL-C)
// 8. replace stencil with foreach for provided iterators
// 9. enclose in namespace
//10. CheckCurrentAutomatizator -> CheckCurrentRowMerger
//10.0 structure/service to handle multiple RowMerger objects in openmp environment
//10.1 user should be able to provide RowMerger when Automatizator is not compiled
//10.2 Automatizator may provide internal structure for RowMerger
#ifdef _MSC_VER
#pragma warning(disable : 4503)
......@@ -1608,5 +1607,5 @@ template<class A> __INLINE INMOST::function_expression<
#endif
#endif
#endif //USE_AUTODIFF
#endif //INMOST_AUTODIFF_ETEXPR_H_INCLUDED
Source/Headers/inmost_mesh.h
View file @ fb24de0b
......@@ -179,7 +179,7 @@ namespace INMOST
public:
reverse_iterator(Mesh * m, const cont_t::reverse_iterator & other ) : cont_t::reverse_iterator(other), m_link(m) {assert(m_link);}
reverse_iterator(const reverse_iterator & other ) : cont_t::reverse_iterator(other), m_link(other.m_link) {assert(m_link);}
ptrdiff_t operator -(const reverse_iterator & other) const {return static_cast<const cont_t::reverse_iterator>(other)-static_cast<const cont_t::reverse_iterator>(*this);}
ptrdiff_t operator -(const reverse_iterator & other) const {return static_cast<const cont_t::reverse_iterator>(*this)-static_cast<const cont_t::reverse_iterator>(other);}
reverse_iterator operator +(size_t n) const {return reverse_iterator(m_link,cont_t::reverse_iterator::operator+(n));}
reverse_iterator operator -(size_t n) const {return reverse_iterator(m_link,cont_t::reverse_iterator::operator-(n));}
reverse_iterator & operator ++() {cont_t::reverse_iterator::operator++(); return *this;}
......@@ -196,6 +196,7 @@ namespace INMOST
public:
const_iterator(Mesh * m, const cont_t::const_iterator & other ) : cont_t::const_iterator(other), m_link(m) {assert(m_link);}
const_iterator(const const_iterator & other ) : cont_t::const_iterator(other), m_link(other.m_link) {assert(m_link);}
ptrdiff_t operator -(const const_iterator & other) const {return static_cast<const cont_t::const_iterator>(*this)-static_cast<const cont_t::const_iterator>(other);}
const_iterator & operator ++() {cont_t::const_iterator::operator++(); return *this;}
const_iterator operator ++(int) {const_iterator ret(*this); cont_t::const_iterator::operator++(); return ret;}
const_iterator & operator --() {cont_t::const_iterator::operator--(); return *this;}
......@@ -210,6 +211,7 @@ namespace INMOST
public:
const_reverse_iterator(Mesh * m, const cont_t::const_reverse_iterator & other) : cont_t::const_reverse_iterator(other), m_link(m) {assert(m_link);}
const_reverse_iterator(const const_reverse_iterator & other ) : cont_t::const_reverse_iterator(other), m_link(other.m_link) {assert(m_link);}
ptrdiff_t operator -(const const_reverse_iterator & other) const {return static_cast<const cont_t::const_reverse_iterator>(*this)-static_cast<const cont_t::const_reverse_iterator>(other);}
const_reverse_iterator & operator ++() {cont_t::const_reverse_iterator::operator++(); return *this;}
const_reverse_iterator operator ++(int) {const_reverse_iterator ret(*this); cont_t::const_reverse_iterator::operator++(); return ret;}
const_reverse_iterator & operator --() {cont_t::const_reverse_iterator::operator--(); return *this;}
......@@ -431,7 +433,7 @@ namespace INMOST
void PrintElementConnectivity() const;
static bool CheckConnectivity (Mesh * m);
//implemented in geometry.cpp
void CastRay (real * pos, real * dir, dynarray< std::pair<Element, real> , 16 > & hits) const;
void CastRay (const real * pos, const real * dir, tiny_map<HandleType, real, 16> & hits) const;
void ComputeGeometricType () const;
void Centroid (real * cnt) const;
void Barycenter (real * cnt) const;
......@@ -604,8 +606,11 @@ namespace INMOST
static Face UniteFaces (ElementArray<Face> & faces, MarkerType del_protect);
static bool TestUniteFaces (const ElementArray<Face> & faces, MarkerType del_protect);
static ElementArray<Face> SplitFace (Face face, const ElementArray<Edge> & edges, MarkerType del_protect); //provide all edges that lay inside face
static bool TestSplitFace (Face face, const ElementArray<Edge> & edges, MarkerType del_protect);
void SwapCells (); //swap back cell and front cell
static bool TestSplitFace (Face face, const ElementArray<Edge> & edges, MarkerType del_protect);
///This function changes Face::BackCell() with Face::FrontCell().
/// \warning Function does not modify normal orientation. You may have
/// to call Face::FixNormalOrientation().
void SwapCells () const;
//implemented in geometry.cpp
Storage::real Area () const;
void Normal (real * nrm) const;
......@@ -1268,8 +1273,8 @@ namespace INMOST
Tag tag_bridge;
Tag tag_redistribute;
private:
void AllocatePrivateMarkers();
void DeallocatePrivateMarkers();
void AllocatePrivateMarkers();
void DeallocatePrivateMarkers();
__INLINE static sparse_rec mkrec (const Tag & t) {sparse_rec ret; ret.tag = t.mem; ret.rec = NULL; return ret;}
__INLINE sparse_type const & MGetSparseLink (integer etypenum, integer ID) const {return GetSparseData(etypenum,links[etypenum][ID]);}
__INLINE sparse_type & MGetSparseLink (integer etypenum, integer ID) {return GetSparseData(etypenum,links[etypenum][ID]);}
......
Source/Headers/inmost_sparse.h
View file @ fb24de0b
......@@ -5,18 +5,44 @@
#include "inmost_common.h"
#if defined(USE_SOLVER)
namespace INMOST
{
namespace Sparse
{
INMOST_MPI_Type GetRowEntryType();
void CreateRowEntryType();
void DestroyRowEntryType();
bool HaveRowEntryType();
#if defined(USE_SOLVER) || defined(USE_AUTODIFF)
/// Retrive MPI type for row entry type
INMOST_MPI_Type GetRowEntryType();
/// Create MPI type for row entry type
void CreateRowEntryType();
/// Release MPI type
void DestroyRowEntryType();
/// Check whether MPI type was created
bool HaveRowEntryType();
/// This class can be used to annotate the matrix
class AnnotationService
{
interval<INMOST_DATA_ENUM_TYPE,std::string> text;
public:
AnnotationService(INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0) { if( end != start ) SetInterval(start,end); }
AnnotationService(const AnnotationService & other) : text(other.text) { }
AnnotationService & operator = (AnnotationService const & other) {text = other.text; return *this;}
~AnnotationService() {}
/// Get the first row index of the distributed matrix interval.
INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return text.get_interval_beg();}
/// Get the last row index of the distributed matrix interval.
INMOST_DATA_ENUM_TYPE GetLastIndex() const {return text.get_interval_end();}
bool Empty() const {return text.empty();}
void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end) { text.set_interval_beg(beg); text.set_interval_end(end); }
void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const {start = text.get_interval_beg(); end = text.get_interval_end();}
std::string & GetAnnotation(INMOST_DATA_ENUM_TYPE row) {assert(!text.empty()); return text[row];}
const std::string & GetAnnotation(INMOST_DATA_ENUM_TYPE row) const {assert(!text.empty()); return text[row];}
void SetAnnotation(INMOST_DATA_ENUM_TYPE row, std::string str) {assert(!text.empty()); text[row] = str;}
};
#endif
#if defined(USE_SOLVER)
/// Distributed vector class.
/// This class can be used to store both local and distributed dense data of real type.
/// For example, to form the right-hand side or initial guess to the solution.
......@@ -85,8 +111,10 @@ namespace INMOST
};
#endif //defined(USE_SOLVER)
#if defined(USE_SOLVER) || defined(USE_AUTODIFF)
/// Class to store the sparse matrix row.
class Row
{
......@@ -162,8 +190,10 @@ namespace INMOST
reverse_iterator rEnd() { return data.rend(); }
const_reverse_iterator rBegin() const { return data.rbegin(); }
const_reverse_iterator rEnd() const { return data.rend(); }
#if defined(USE_SOLVER)
/// Return the scalar product of the current sparse row by a dense Vector.
INMOST_DATA_REAL_TYPE RowVec(Vector & x) const; // returns A(row) * x
#endif
void MoveRow(Row & new_pos) {data = new_pos.data;} //here move constructor and std::move may be used in future
/// Set the vector entries by zeroes.
void Zero() {for(iterator it = Begin(); it != End(); ++it) it->second = 0;}
......@@ -185,7 +215,9 @@ namespace INMOST
static void MergeSortedRows(INMOST_DATA_REAL_TYPE alpha, const Row & left, INMOST_DATA_REAL_TYPE beta, const Row & right, Row & output);
};
#endif //defined(USE_SOLVER) || defined(USE_AUTODIFF)
#if defined(USE_SOLVER) || defined(USE_AUTODIFF)
/// Class to store the compressed symmetric matrix of a hessian row.
class HessianRow
{
......@@ -298,7 +330,11 @@ namespace INMOST
static void MergeJacobianHessian(INMOST_DATA_REAL_TYPE a, const Row & JL, const Row & JR, INMOST_DATA_REAL_TYPE b, const HessianRow & HL, INMOST_DATA_REAL_TYPE c, const HessianRow & HR, HessianRow & output);
static void MergeJacobianHessian(INMOST_DATA_REAL_TYPE a, const Row & JL, const Row & JR, INMOST_DATA_REAL_TYPE b, const HessianRow & H, HessianRow & output);
};
#endif //defined(USE_SOLVER) || defined(USE_AUTODIFF)
#if defined(USE_SOLVER)
/// This class can be used for shared access to matrix with OpenMP.
class LockService
{
......@@ -322,29 +358,6 @@ namespace INMOST
bool UnLock(INMOST_DATA_ENUM_TYPE row);
};
/// This class can be used to annotate the matrix
class AnnotationService
{
interval<INMOST_DATA_ENUM_TYPE,std::string> text;
public:
AnnotationService(INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0) { if( end != start ) SetInterval(start,end); }
AnnotationService(const AnnotationService & other) : text(other.text) { }
AnnotationService & operator = (AnnotationService const & other) {text = other.text; return *this;}
~AnnotationService() {}
/// Get the first row index of the distributed matrix interval.
INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return text.get_interval_beg();}
/// Get the last row index of the distributed matrix interval.
INMOST_DATA_ENUM_TYPE GetLastIndex() const {return text.get_interval_end();}
bool Empty() const {return text.empty();}
void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end) { text.set_interval_beg(beg); text.set_interval_end(end); }
void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const {start = text.get_interval_beg(); end = text.get_interval_end();}
std::string & GetAnnotation(INMOST_DATA_ENUM_TYPE row) {assert(!text.empty()); return text[row];}
const std::string & GetAnnotation(INMOST_DATA_ENUM_TYPE row) const {assert(!text.empty()); return text[row];}
void SetAnnotation(INMOST_DATA_ENUM_TYPE row, std::string str) {assert(!text.empty()); text[row] = str;}
};
/// Class to store the distributed sparse matrix by compressed rows.
/// The format used to store sparse matrix is analogous to Compressed Row Storage format (CRS).
/// @see http://netlib.org/linalg/html_templates/node91.html
......@@ -414,7 +427,11 @@ namespace INMOST
/// Get the matrix name specified in the main constructor.
std::string GetName() const {return name;}
};
#endif //defined(USE_SOLVER)
#if defined(USE_SOLVER)
/// Class to store the distributed sparse hessian hyper matrix by compressed symmetric matrices.
class HessianMatrix
{
......@@ -482,7 +499,11 @@ namespace INMOST
/// Get the matrix name specified in the main constructor.
std::string GetName() const {return name;}
};
#endif //defined(USE_SOLVER)
#if defined(USE_SOLVER) || defined(USE_AUTODIFF)
/// This class may be used to sum multiple sparse rows.
/// \warning
/// In parallel column indices of the matrix may span wider then
......@@ -538,10 +559,12 @@ namespace INMOST
/// @param interval_end Last index in linked list.
/// @param Sorted Result should be sorted or not.
RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true);
#if defined(USE_SOLVER)
/// Constructor that gets sizes from the matrix.
/// @param A Matrix to get sizes from.
/// @param Sorted Result should be sorted.
RowMerger(Matrix & A, bool Sorted = true);
#endif //USE_SOLVER
/// Destructor.
~RowMerger();
/// Resize linked list for new interval.
......@@ -552,6 +575,7 @@ namespace INMOST
/// @param interval_end Last index in linked list.
/// @param Sorted Result should be sorted or not.
void Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted = true);
#if defined(USE_SOLVER)
/// Resize linked list for new matrix.
/// \warning
/// All contents of linked list will be lost after resize.
......@@ -559,6 +583,7 @@ namespace INMOST
/// @param A Matrix to get sizes from.
/// @param Sorted Result should be sorted or not.
void Resize(Matrix & A, bool Sorted = true);
#endif //USE_SOLVER
/// Clear linked list.
void Clear();
/// Add a row with a coefficient into empty linked list.
......@@ -614,8 +639,8 @@ namespace INMOST
iterator Begin() {return iterator(&LinkedList);}
iterator End() {iterator ret(&LinkedList); ret.pos = EOL; return ret;}
};
}
}
#endif
#endif //defined(USE_SOLVER) || defined(USE_AUTODIFF)
} //namespace Sparse
} //namespace INMOST
#endif
\ No newline at end of file
#endif //INMOST_SPARSE_INCLUDED
\ No newline at end of file
Source/Headers/inmost_variable.h
View file @ fb24de0b
......@@ -9,11 +9,7 @@
#include <sstream> //for debug
#include <new>
#if defined(USE_AUTODIFF) && (!defined(USE_MESH) || !defined(USE_SOLVER))
#warning "USE_AUTODIFF require USE_MESH"
#undef USE_AUTODIFF
#endif
#if defined(USE_AUTODIFF) && defined(USE_MESH)
//TODO:
// 1. Incorporate tables
......@@ -32,7 +28,6 @@
#if defined(USE_AUTODIFF)
namespace INMOST
{
......@@ -582,9 +577,10 @@ template<class A> __INLINE
return INMOST::stencil_expression<typename A::Var>(tmp);
}
#endif //defined(USE_AUTODIFF) && defined(USE_MESH)
#endif
#endif
#endif //INMOST_AUTODIFF_ETVAR_H_INCLUDED
Source/IO/mesh_ecl_file.cpp
View file @ fb24de0b
This source diff could not be displayed because it is too large. You can view the blob instead.
Source/Mesh/geometry.cpp
View file @ fb24de0b
This diff is collapsed.
Source/Mesh/incident_matrix.hpp
View file @ fb24de0b
......@@ -51,6 +51,38 @@ namespace INMOST
}
} orient_face;
class AbstractCoords
{
public:
virtual Storage::real_array Get(const Node & s) const = 0;
virtual AbstractCoords * Copy() const = 0;
virtual ~AbstractCoords() {}
};
class GridCoords : public AbstractCoords
{
public:
Storage::real_array Get(const Node & s) const
{
return s.Coords();
}
AbstractCoords * Copy() const {return new GridCoords(*this);}
};
class TagCoords : public AbstractCoords
{
Tag t;
public:
TagCoords(Tag _t) : t(_t) {}
TagCoords & operator =(TagCoords const & b) {t = b.t; return *this;}
TagCoords(const TagCoords & b) : t(b.t) {}
Storage::real_array Get(const Node & s) const
{
return s.RealArray(t);
}
AbstractCoords * Copy() const {return new TagCoords(*this);}
};
template<class T>
class incident_matrix
{
......@@ -70,6 +102,7 @@ namespace INMOST
//std::vector< std::pair<std::vector<int>,double> > remember;
double min_loop_measure;
bool print;
AbstractCoords * coords;
bool do_hide_row(unsigned k)
{
......@@ -130,7 +163,7 @@ namespace INMOST
if( data[0]->GetElementDimension() == 1 ) //this is edge //use geometric dimension here for 2d compatibility
{
//calculate area
int mdim = data[0]->GetMeshLink()->GetDimensions();
//int mdim = data[0]->GetMeshLink()->GetDimensions();
ElementArray<Node> nodes,n1,n2;
n1 = data[0]->getNodes();
n2 = data[1]->getNodes();
......@@ -155,22 +188,12 @@ namespace INMOST
}
Storage::real x[3] = {0,0,0};
Storage::real x0[3] = {0,0,0};
for(unsigned i = 0; i < nodes.size(); i++)
{
Storage::real_array v = nodes[i].Coords();
x0[0] += v[0];
x0[1] += v[1];
if( mdim == 3 ) x0[2] += v[2];
}
x0[0] /= (Storage::real)nodes.size();
x0[1] /= (Storage::real)nodes.size();
x0[2] /= (Storage::real)nodes.size();
for(unsigned i = 0; i < nodes.size()-1; i++)
Storage::real_array x0 = coords->Get(nodes[0]);
for(unsigned i = 1; i < nodes.size()-1; i++)
{
Storage::real_array v1 = nodes[i].Coords();
Storage::real_array v2 = nodes[i+1].Coords();
if( mdim == 3 )
Storage::real_array v1 = coords->Get(nodes[i]);
Storage::real_array v2 = coords->Get(nodes[i+1]);
if( v1.size() == 3 && v2.size() == 3 )
{
x[0] += (v1[1]-x0[1])*(v2[2]-x0[2]) - (v1[2]-x0[2])*(v2[1]-x0[1]);
x[1] += (v1[2]-x0[2])*(v2[0]-x0[0]) - (v1[0]-x0[0])*(v2[2]-x0[2]);
......@@ -271,11 +294,34 @@ namespace INMOST
} while(true);
data.RemPrivateMarker(mrk);
mesh->ReleasePrivateMarker(mrk);
Storage::real cnt[3], nrm[3];
Storage::real nrm[3], cnt[3];
Storage::real_array v0,v1,v2;
for(typename ElementArray<T>::size_type j = 0; j < data.size(); j++)
{
data[j]->Centroid(cnt);
data[j]->getAsFace()->Normal(nrm);
//compute normal to face
ElementArray<Node> nodes = data[j]->getAsFace()->getNodes();
nrm[0] = nrm[1] = nrm[2] = 0;
v0 = v1 = coords->Get(nodes[0]);
cnt[0] = v1[0];
cnt[1] = v1[1];
cnt[2] = v1[2];
for(unsigned i = 0; i < nodes.size(); i++)
{
v2 = coords->Get(nodes[(i+1)%nodes.size()]);
cnt[0] += v2[0];
cnt[1] += v2[1];
cnt[2] += v2[2];
nrm[0] += (v1[1]-v0[1])*(v2[2]-v0[2]) - (v1[2]-v0[2])*(v2[1]-v0[1]);
nrm[1] += (v1[2]-v0[2])*(v2[0]-v0[0]) - (v1[0]-v0[0])*(v2[2]-v0[2]);
nrm[2] += (v1[0]-v0[0])*(v2[1]-v0[1]) - (v1[1]-v0[1])*(v2[0]-v0[0]);
v1.swap(v2);
}
nrm[0] *= 0.5;
nrm[1] *= 0.5;
nrm[2] *= 0.5;
cnt[0] /= (Storage::real)nodes.size();
cnt[1] /= (Storage::real)nodes.size();
cnt[2] /= (Storage::real)nodes.size();
measure += (data[j]->GetPrivateMarker(rev) ? -1.0 : 1.0)*dot_prod(cnt,nrm);
}
data.RemPrivateMarker(rev);
......@@ -428,8 +474,8 @@ namespace INMOST
}
}
template<typename InputIterator>
incident_matrix(Mesh * mesh, InputIterator beg, InputIterator end, typename ElementArray<T>::size_type num_inner, bool print = false)
: mesh(mesh), head_column(beg,end), min_loop(), print(print)
incident_matrix(Mesh * mesh, InputIterator beg, InputIterator end, typename ElementArray<T>::size_type num_inner, const AbstractCoords & _coords = GridCoords(), bool print = false)
: mesh(mesh), head_column(beg,end), min_loop(), print(print), coords(_coords.Copy())
{
min_loop.SetMeshLink(mesh);
temp_loop.SetMeshLink(mesh);
......@@ -494,7 +540,7 @@ namespace INMOST
: mesh(other.mesh), matrix(other.matrix), head_column(other.head_column), head_row(other.head_row),
head_row_count(other.head_row_count), min_loop(other.min_loop),
hide_row(other.hide_row), hide_column(other.hide_column),
stub_row(other.stub_row)
stub_row(other.stub_row), coords(other.coords.Copy())
{
}
incident_matrix & operator =(const incident_matrix & other)
......@@ -508,10 +554,12 @@ namespace INMOST
hide_row = other.hide_row;
hide_column = other.hide_column;
stub_row = other.stub_row;
coords = other.coords.Copy();
return *this;
}
~incident_matrix()
{
delete coords;
}
bool find_shortest_loop(ElementArray<T> & ret)
{
......
Source/Mesh/mesh.cpp
View file @ fb24de0b
......@@ -1309,52 +1309,37 @@ namespace INMOST
*/
case Element::Hex:
{
MarkerType mrk = CreatePrivateMarker();
MarkerType cemrk = CreatePrivateMarker();
MarkerType femrk = CreatePrivateMarker();
//printf("%lx %lx %lx\n",mrk,cemrk,femrk);
ElementArray<Face> faces = c->getFaces();
Face face = faces[0];
ret.reserve(8);
ElementArray<Node> verts = face->getNodes();
if( face->BackCell() == c )
for(ElementArray<Node>::iterator it = verts.begin(); it != verts.end(); it++)
{
ret.push_back(*it);
it->SetPrivateMarker(mrk);
}
ret.insert(ret.end(),verts.begin(),verts.end());
else
for(ElementArray<Node>::reverse_iterator it = verts.rbegin(); it != verts.rend(); it++)
{
ret.push_back(*it);
it->SetPrivateMarker(mrk);
}
ret.insert(ret.end(),verts.rbegin(),verts.rend());
ElementArray<Edge> c_edges = c->getEdges();
for(ElementArray<Edge>::iterator it = c_edges.begin(); it != c_edges.end(); it++)
it->SetPrivateMarker(cemrk);
ElementArray<Edge> f_edges = face->getEdges();
for(ElementArray<Edge>::iterator it = f_edges.begin(); it != f_edges.end(); it++)
it->SetPrivateMarker(femrk);
c_edges.SetPrivateMarker(cemrk);
f_edges.SetPrivateMarker(femrk);
for(unsigned int k = 0; k < 4; k++)
{
ElementArray<Edge> v_edges = ret[k]->getEdges();
ElementArray<Edge> v_edges = ret[k]->getEdges(cemrk);
for(ElementArray<Edge>::iterator it = v_edges.begin(); it != v_edges.end(); it++)
{
if( it->GetPrivateMarker(cemrk) && !it->GetPrivateMarker(femrk) )
if( !it->GetPrivateMarker(femrk) )
{
ElementArray<Node> nn = it->getNodes();
if( nn[0].GetPrivateMarker(mrk) )
ret.push_back(nn[1]);
if( it->getBeg() == ret[k] )
ret.push_back(it->getEnd());
else
ret.push_back(nn[0]);
ret.push_back(it->getBeg());
break;
}
}
ret[k]->RemPrivateMarker(mrk);
}
for(ElementArray<Edge>::iterator it = c_edges.begin(); it != c_edges.end(); it++) it->RemPrivateMarker(cemrk);
for(ElementArray<Edge>::iterator it = f_edges.begin(); it != f_edges.end(); it++) it->RemPrivateMarker(femrk);
ReleasePrivateMarker(mrk);
c_edges.RemPrivateMarker(cemrk);
f_edges.RemPrivateMarker(femrk);
ReleasePrivateMarker(cemrk);
ReleasePrivateMarker(femrk);
break;
......@@ -1367,7 +1352,6 @@ namespace INMOST
*/
case Element::Prism:
{
MarkerType mrk = CreatePrivateMarker();
MarkerType cemrk = CreatePrivateMarker();
MarkerType femrk = CreatePrivateMarker();
ret.reserve(6);
......@@ -1381,39 +1365,28 @@ namespace INMOST
}
ElementArray<Node> verts = face->getNodes();
if( face->BackCell() == c )
for(ElementArray<Node>::iterator it = verts.begin(); it != verts.end(); it++)
{
ret.push_back(*it);
it->SetPrivateMarker(mrk);
}
ret.insert(ret.end(),verts.begin(),verts.end());
else
for(ElementArray<Node>::reverse_iterator it = verts.rbegin(); it != verts.rend(); it++)
{
ret.push_back(*it);
it->SetPrivateMarker(mrk);
}
ret.insert(ret.end(),verts.rbegin(),verts.rend());
ElementArray<Edge> c_edges = c->getEdges();
for(ElementArray<Edge>::iterator it = c_edges.begin(); it != c_edges.end(); it++) it->SetPrivateMarker(cemrk);
ElementArray<Edge> f_edges = face->getEdges();
for(ElementArray<Edge>::iterator it = f_edges.begin(); it != f_edges.end(); it++) it->SetPrivateMarker(femrk);
c_edges.SetPrivateMarker(cemrk);
f_edges.SetPrivateMarker(femrk);
for(ElementArray<Cell>::size_type k = 0; k < 3; k++)
{
ElementArray<Edge> v_edges = ret[k]->getEdges();
ElementArray<Edge> v_edges = ret[k]->getEdges(cemrk);
for(ElementArray<Edge>::iterator it = v_edges.begin(); it != v_edges.end(); it++)
if( it->GetPrivateMarker(cemrk) && !it->GetPrivateMarker(femrk) )
if( !it->GetPrivateMarker(femrk) )
{
ElementArray<Node> nn = it->getNodes();
if( nn[0].GetPrivateMarker(mrk) )
ret.push_back(nn[1]);
if( it->getBeg() == ret[k] )
ret.push_back(it->getEnd());
else
ret.push_back(nn[0]);
ret.push_back(it->getBeg());
break;
}
ret[k]->RemPrivateMarker(mrk);
}
for(ElementArray<Edge>::iterator it = c_edges.begin(); it != c_edges.end(); it++) it->RemPrivateMarker(cemrk);
for(ElementArray<Edge>::iterator it = f_edges.begin(); it != f_edges.end(); it++) it->RemPrivateMarker(femrk);
ReleasePrivateMarker(mrk);
c_edges.RemPrivateMarker(cemrk);
f_edges.RemPrivateMarker(femrk);
ReleasePrivateMarker(cemrk);
ReleasePrivateMarker(femrk);
break;
......@@ -1527,13 +1500,12 @@ namespace INMOST
SetPrivateMarker(ilc[j],mrk);
}
}
for(Element::adj_type::size_type i = 0; i < lc.size(); i++) if( !Hidden(lc[i]) ) //iterate over faces
{
Element::adj_type & ilc = LowConn(lc[i]);
for(Element::adj_type::size_type j = 0; j < ilc.size(); j++) if( !Hidden(ilc[j]) ) //iterate over face edges
RemPrivateMarker(ilc[j],mrk);