diff --git a/Examples/GridTools/CMakeLists.txt b/Examples/GridTools/CMakeLists.txt index 4fb29531bc2e39d3f313c246322098259a8a1efc..633dcd893aafd6fa09320022fa3ad84d59d0aa1f 100644 --- a/Examples/GridTools/CMakeLists.txt +++ b/Examples/GridTools/CMakeLists.txt @@ -6,6 +6,7 @@ add_executable(UniteFaces unite_faces.cpp) add_executable(Dual dual.cpp) add_executable(Slice slice.cpp) add_executable(SliceFunc slice_func.cpp) +add_executable(SplitNonplanar split_nonplanar.cpp) target_link_libraries(FixFaults inmost) if(USE_MPI) @@ -71,6 +72,15 @@ endif(USE_MPI) install(TARGETS SliceFunc EXPORT inmost-targets RUNTIME DESTINATION bin) +target_link_libraries(SplitNonplanar inmost) +if(USE_MPI) + message("linking SplitNonplanar with MPI") + target_link_libraries(SplitNonplanar ${MPI_LIBRARIES}) + if(MPI_LINK_FLAGS) + set_target_properties(SplitNonplanar PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") + endif() +endif(USE_MPI) +install(TARGETS SplitNonplanar EXPORT inmost-targets RUNTIME DESTINATION bin) diff --git a/Examples/GridTools/dual.cpp b/Examples/GridTools/dual.cpp index 1c3bc0e6652f3902267030536a1acc42b878a801..382aed0c95154a30a0bde9f52e758fecaa2994c8 100644 --- a/Examples/GridTools/dual.cpp +++ b/Examples/GridTools/dual.cpp @@ -19,7 +19,7 @@ int main(int argc, char ** argv) Mesh::GeomParam table; table[ORIENTATION] = FACE; - table[CENTROID] = FACE | CELL; + table[BARYCENTER] = FACE | CELL | EDGE; A.PrepareGeometricData(table); @@ -28,13 +28,13 @@ int main(int argc, char ** argv) for(Mesh::iteratorCell it = A.BeginCell(); it != A.EndCell(); ++it) { real cnt[3]; - it->Centroid(cnt); + it->Barycenter(cnt); it->RemoteReference(cell2node) = RemoteHandleType(&B,B.CreateNode(cnt).GetHandle()); } for(Mesh::iteratorFace it = A.BeginFace(); it != A.EndFace(); ++it) if( it->Boundary() ) { real cnt[3]; - it->Centroid(cnt); + it->Barycenter(cnt); it->RemoteReference(cell2node) = RemoteHandleType(&B,B.CreateNode(cnt).GetHandle()); } //add corner nodes @@ -54,7 +54,7 @@ int main(int argc, char ** argv) { corners++; it->SetMarker(corner); - it->Centroid(cnt); + it->Barycenter(cnt); it->RemoteReference(cell2node) = RemoteHandleType(&B,B.CreateNode(cnt).GetHandle()); } } @@ -66,7 +66,7 @@ int main(int argc, char ** argv) { corners++; it->SetMarker(corner); - it->Centroid(cnt); + it->Barycenter(cnt); it->RemoteReference(cell2node) = RemoteHandleType(&B,B.CreateNode(cnt).GetHandle()); } else if( ncorners == 2 ) @@ -75,9 +75,9 @@ int main(int argc, char ** argv) Node n1 = edges[0]->getBeg() == it->self() ? edges[0]->getEnd() : edges[0]->getBeg(); Node n2 = edges[1]->getBeg() == it->self() ? edges[1]->getEnd() : edges[1]->getBeg(); real cnt1[3], cnt2[3], r1[3],r2[3], l, dot; - it->Centroid(cnt); - n1->Centroid(cnt1); - n2->Centroid(cnt2); + it->Barycenter(cnt); + n1->Barycenter(cnt1); + n2->Barycenter(cnt2); r1[0] = cnt[0] - cnt1[0]; r1[1] = cnt[1] - cnt1[1]; r1[2] = cnt[2] - cnt1[2]; @@ -98,7 +98,7 @@ int main(int argc, char ** argv) { corners++; it->SetMarker(corner); - it->Centroid(cnt); + it->Barycenter(cnt); it->RemoteReference(cell2node) = RemoteHandleType(&B,B.CreateNode(cnt).GetHandle()); } } diff --git a/Examples/GridTools/split_nonplanar.cpp b/Examples/GridTools/split_nonplanar.cpp new file mode 100644 index 0000000000000000000000000000000000000000..eb2da34d16737b3940f330cf9dab74b81df4e47a --- /dev/null +++ b/Examples/GridTools/split_nonplanar.cpp @@ -0,0 +1,380 @@ +#include +#include "inmost.h" + + +using namespace INMOST; +typedef Storage::real real; +typedef Storage::real_array real_array; + +real dot(real a[3], real b[3]) +{ + return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; +} + +void cross(real a[3], real b[3], real out[3]) +{ + out[0] = a[1]*b[2] - a[2]*b[1]; + out[1] = a[2]*b[0] - a[0]*b[2]; + out[2] = a[0]*b[1] - a[1]*b[0]; +} + +static void normalize(double v[3]) +{ + double d = sqrt(dot(v, v)); + if (d) for (int k = 0; k < 3; ++k) v[k] /= d; +} + + +struct node +{ + real coo[2]; + node(const real _coo[2]) + { + coo[0] = _coo[0]; + coo[1] = _coo[1]; + } + node(const node & b) + { + coo[0] = b.coo[0]; + coo[1] = b.coo[1]; + } + node & operator = (node const & b) + { + coo[0] = b.coo[0]; + coo[1] = b.coo[1]; + return *this; + } +}; + + +struct segment +{ + real beg[2]; + real end[2]; + segment(const real _beg[2], const real _end[2]) + { + for(int k = 0; k < 2; ++k) + { + beg[k] = _beg[k]; + end[k] = _end[k]; + } + } + segment(const segment & b) + { + for(int k = 0; k < 2; ++k) + { + beg[k] = b.beg[k]; + end[k] = b.end[k]; + } + } + segment & operator = (segment const & b) + { + for(int k = 0; k < 2; ++k) + { + beg[k] = b.beg[k]; + end[k] = b.end[k]; + } + return *this; + } +}; + + +//check whether two segments intersect +bool intersect_segments(const segment & a, const segment & b, bool print = false) +{ + const real eps = 1.0e-9; + real div1 = (a.beg[0] - a.end[0])*(b.beg[1] - b.end[1]); + real div2 = (a.beg[1] - a.end[1])*(b.beg[0] - b.end[0]); + real div = div1 - div2; + real t, find[2]; + if( print ) + { + std::cout << "(" << a.beg[0] << "," << a.beg[1] << ",0) <-> (" << a.end[0] << "," << a.end[1] << ",0)" << std::endl; + std::cout << "(" << b.beg[0] << "," << b.beg[1] << ",0) <-> (" << b.end[0] << "," << b.end[1] << ",0)" << std::endl; + } + if( print ) std::cout << "div = " << div << " div1 " << div1 << " div2 " << div2 << " against " << eps*(std::abs(div1) + std::abs(div2)) << std::endl; + if (std::abs(div) <= eps*(std::abs(div1) + std::abs(div2))) //segments are parallel + { + //return true; + if( print ) std::cout << "parallel segments, div = " << div << std::endl; + if( std::abs(a.beg[0]-a.end[0]) < eps*(std::abs(a.beg[0])+std::abs(a.end[0])) ) //parallel to x + { + real a0 = a.beg[1], a1 = a.end[1]; + real b0 = b.beg[1], b1 = b.end[1]; + if( print ) std::cout << "parallel to x, projection on Oy: a [" << a0 << "," << a1 << "], b [" << b0 << "," << b1 << "]" << std::endl; + if( a0 > a1 ) std::swap(a0,a1); + if( b0 > b1 ) std::swap(b0,b1); + if( a1 > b1 ) //change a and b + { + std::swap(a0,b0); + std::swap(a1,b1); + } + if( a1 > b0+eps ) + { + if( print ) std::cout << "intersection" << std::endl; + return true; + } + if( print ) std::cout << "no intersection" << std::endl; + return false; + } + //reconstruct k in y = kx+b + real ak = (a.end[1]-a.beg[1])/(a.end[0]-a.beg[0]); + real bk = (b.end[1]-b.beg[1])/(b.end[0]-b.beg[0]); + //reconstruct bs in y = kx+b + real ab = (a.end[0]*a.beg[1] - a.beg[0]*a.end[1])/(a.end[0]-a.beg[0]); + real bb = (b.end[0]*b.beg[1] - b.beg[0]*b.end[1])/(b.end[0]-b.beg[0]); + if( print ) + { + std::cout << "parallel lines:" << std::endl; + std::cout << "a: y = " << ak << "*x + " << ab << std::endl; + std::cout << "b: y = " << bk << "*x + " << bb << std::endl; + std::cout << "kdiff: " << ak - bk << std::endl; + std::cout << "bdiff: " << ab - bb << " against " << eps*(std::abs(ab)+std::abs(bb)) << std::endl; + } + if( std::abs(ab-bb) <= eps*(std::abs(ab)+std::abs(bb)) ) + { + //test segments + int ind = ak > 1.0 ? 1 : 0; //choose projection of segments + real a0 = a.beg[ind], a1 = a.end[ind]; + real b0 = b.beg[ind], b1 = b.end[ind]; + if( print ) std::cout << "projection on O" << (ind?"y":"x") << ": a [" << a0 << "," << a1 << "], b [" << b0 << "," << b1 << "]" << std::endl; + if( a0 > a1 ) std::swap(a0,a1); + if( b0 > b1 ) std::swap(b0,b1); + if( a1 > b1 ) //change a and b + { + std::swap(a0,b0); + std::swap(a1,b1); + } + if( a1 > b0+eps ) + { + if( print ) std::cout << "intersection" << std::endl; + return true; + } + if( print ) std::cout << "no intersection" << std::endl; + return false; + } + //printf("divisor is zero\n"); + return false; + } + find[0] = ((a.beg[0]*a.end[1] - a.beg[1]*a.end[0])*(b.beg[0] - b.end[0]) - (a.beg[0] - a.end[0])*(b.beg[0]*b.end[1] - b.beg[1]*b.end[0])) / div; + find[1] = ((a.beg[0]*a.end[1] - a.beg[1]*a.end[0])*(b.beg[1] - b.end[1]) - (a.beg[1] - a.end[1])*(b.beg[0]*b.end[1] - b.beg[1]*b.end[0])) / div; + if( print ) std::cout << "intersection: (" << find[0] << "," << find[1] << ")" << std::endl; + //probably some of these tests are redundant + bool inter = true; + if (std::abs(a.end[0] - a.beg[0]) > eps) + { + t = (find[0] - a.beg[0]) / (a.end[0] - a.beg[0]); + if( print ) std::cout << "coef a on Ox " << t << std::endl; + if (t < eps || t > 1.0 - eps) inter = false; + } + if (std::abs(a.end[1] - a.beg[1]) > eps) + { + t = (find[1] - a.beg[1]) / (a.end[1] - a.beg[1]); + if( print ) std::cout << "coef a on Oy " << t << std::endl; + if (t < eps || t > 1.0 - eps) inter = false; + } + if (std::abs(b.end[0] - b.beg[0]) > eps) + { + t = (find[0] - b.beg[0]) / (b.end[0] - b.beg[0]); + if( print ) std::cout << "coef b on Ox " << t << std::endl; + if (t < eps || t > 1.0 - eps) inter = false; + } + if (std::abs(b.end[1] - b.beg[1]) > eps) + { + t = (find[1] - b.beg[1]) / (b.end[1] - b.beg[1]); + if( print ) std::cout << "coef b on Oy " << t << std::endl; + if (t < eps || t > 1.0 - eps) inter = false; + } + if( inter ) + { + if( print ) std::cout << "intersection" << std::endl; + return true; + } + if( print ) std::cout << "no intersection" << std::endl; + return false; +} + +//check whether any of the segments from setb intersect any segments from seta +bool check_intersect(const std::vector & segmentsa, const std::vector & segmentsb, bool print = false) +{ + if( print ) + { + std::cout << "segments [" << segmentsa.size() << "]:" << std::endl; + for(int k = 0; k < segmentsa.size(); ++k) + std::cout << "(" << segmentsa[k].beg[0] << "," << segmentsa[k].beg[1] << ",0) <-> (" << segmentsa[k].end[0] << "," << segmentsa[k].end[1] << ",0)" << std::endl; + std::cout << "joints [" << segmentsb.size() << "]:" << std::endl; + for(int k = 0; k < segmentsb.size(); ++k) + std::cout << "(" << segmentsb[k].beg[0] << "," << segmentsb[k].beg[1] << ",0) <-> (" << segmentsb[k].end[0] << "," << segmentsb[k].end[1] << ",0)" << std::endl; + } + for(int i = 0; i < (int)segmentsb.size(); ++i) + { + for(int j = 0; j < (int)segmentsa.size(); ++j) + { + if( print ) std::cout << "check segment " << j << " and joint " << i << std::endl; + if( intersect_segments(segmentsa[j],segmentsb[i],print) ) + { + if( print ) std::cout << "there is intersection of segment " << j << " and joint " << i << std::endl; + return true; + } + } + } + if( print ) std::cout << "there is no intersection" << std::endl; + return false; +} + + +//check whether any of the centers of setb drop outside of loop represented by seta +bool check_dropout(const std::vector & loop, const std::vector & segmentsb) +{ + for(int i = 0; i < (int)segmentsb.size(); ++i) + { + real cnt[2]; + cnt[0] = (segmentsb[i].beg[0]+segmentsb[i].end[0])*0.5; + cnt[1] = (segmentsb[i].beg[1]+segmentsb[i].end[1])*0.5; + + int cn = 0; // the crossing number counter + // loop through all edges of the polygon + for (int j = 0; j < loop.size()-1; j++) // edge from V[i] to V[i+1] + { + real px = cnt[0]; + real py = cnt[1]; + real v0x = loop[j].coo[0]; + real v0y = loop[j].coo[1]; + real v1x = loop[j+1].coo[0]; + real v1y = loop[j+1].coo[1]; + if ( ((v0y <= py) && (v1y > py)) // an upward crossing + || ((v0y > py) && (v1y <= py)))// a downward crossing + { + // compute the actual edge-ray intersect x-coordinate + real vt = (py - v0y) / (v1y - v0y); + if (px < v0x + vt * (v1x - v0x)) // P.x < intersect + ++cn; // a valid crossing of y=P.y right of P.x + } + } + if( !(cn&1) ) //out of polygon + return true; + } + return false; +} + + +int main(int argc, char ** argv) +{ + if( argc < 2 ) + { + printf("Usage: %s input_mesh [output_mesh]\n",argv[0]); + return -1; + } + + Mesh A(""); + A.Load(argv[1]); + + Mesh::GeomParam table; + table[ORIENTATION] = FACE; + table[BARYCENTER] = FACE | CELL | EDGE; + A.PrepareGeometricData(table); + + Tag proj_coords = A.CreateTag("PROJECTED_COORDS",DATA_REAL,NODE,NONE,2); + + real nrm[3], cnt[3], ncnt[3], scnt[3], pcnt[3], fcnt[3], ray[3], orthx[3], orthy[3], d, nd, c[2]; + for(Mesh::iteratorFace it = A.BeginFace(); it != A.EndFace(); ++it) + { + if( !it->Planarity() ) + { + std::vector segments; + it->UnitNormal(nrm); + it->Centroid(cnt); + d = dot(cnt,nrm); + ElementArray face_nodes = it->getNodes(); + //project coordinates + for(ElementArray::iterator mt = face_nodes.begin(); mt != face_nodes.end(); ++mt) + { + mt->Centroid(ncnt); + nd = dot(ncnt,nrm)-d; + //project onto plane + pcnt[0] = ncnt[0] - cnt[0]; + pcnt[1] = ncnt[1] - cnt[1]; + pcnt[2] = ncnt[2] - cnt[2]; + //choose the basis + if( mt == face_nodes.begin() ) + { + orthx[0] = pcnt[0]; + orthx[1] = pcnt[1]; + orthx[2] = pcnt[2]; + normalize(orthx); + cross(orthx,nrm,orthy); + } + //compute and record the coords + real_array c = mt->RealArray(proj_coords); + c[0] = dot(orthx,pcnt); + c[1] = dot(orthy,pcnt); + } + + //gather all the segments in projected coordinates + ElementArray face_edges = it->getEdges(); + for(ElementArray::iterator mt = face_edges.begin(); mt != face_edges.end(); ++mt) + segments.push_back(segment(mt->getBeg().RealArray(proj_coords).data(),mt->getEnd().RealArray(proj_coords).data())); + //form a loop out of face nodes + std::vector loop; + for(ElementArray::iterator mt = face_nodes.begin(); mt != face_nodes.end(); ++mt) + loop.push_back(node(mt->RealArray(proj_coords).data())); + loop.push_back(loop.front()); + //try each adjacent node + bool success = false; + for(ElementArray::iterator jt = face_nodes.begin(); jt != face_nodes.end() && !success; ++jt) + { + std::vector joints; + //check that the node sees centers of all segments + for(ElementArray::iterator mt = face_edges.begin(); mt != face_edges.end(); ++mt) + { + if( mt->getBeg() != jt->self() && mt->getEnd() != jt->self() ) //skip adjacent to the node + { + scnt[0] = (mt->getBeg()->RealArray(proj_coords)[0]+mt->getEnd()->RealArray(proj_coords)[0])*0.5; + scnt[1] = (mt->getBeg()->RealArray(proj_coords)[1]+mt->getEnd()->RealArray(proj_coords)[1])*0.5; + joints.push_back(segment(jt->RealArray(proj_coords).data(),scnt)); + } + } + + if( !check_intersect(segments,joints) && !check_dropout(loop,joints) ) + { + //found a suitable node + ElementArray new_edges(&A); + ElementArray edge_nodes(&A,2); + edge_nodes[0] = jt->self(); + for(ElementArray::iterator kt = face_nodes.begin(); kt != face_nodes.end(); ++kt) + { + if( kt->self() != jt->self() ) + { + edge_nodes[1] = kt->self(); + std::pair edge = A.CreateEdge(edge_nodes); + if( edge.second ) //this is actually a new edge + new_edges.push_back(edge.first); + } + } + ElementArray split_faces = Face::SplitFace(it->self(),new_edges,0); + success = true; + } + } + } + } + + if( A.HaveTag("GRIDNAME") ) + { + Storage::bulk_array nameA = A.self().BulkArray(A.GetTag("GRIDNAME")); + std::string ins = "_split_nonplanar"; + nameA.insert(nameA.end(),ins.c_str(),ins.c_str()+ins.size()); + } + + if( argc > 2 ) + { + std::cout << "Save to " << argv[2] << std::endl; + A.Save(argv[2]); + } + else + { + std::cout << "Save to out.vtk" << std::endl; + A.Save("out.vtk"); + } + + return 0; +} \ No newline at end of file diff --git a/Source/Headers/inmost_dense.h b/Source/Headers/inmost_dense.h index 9afc81fc17ab8da5f73c0f5eaf9e01f3e1f3a8b0..0de59c1370036a0007d14d6987464c18910febe9 100644 --- a/Source/Headers/inmost_dense.h +++ b/Source/Headers/inmost_dense.h @@ -1483,6 +1483,7 @@ namespace INMOST return *this; } + template template Matrix::type> @@ -1503,6 +1504,43 @@ namespace INMOST return ret; } + + /* + template + template + Matrix::type> + AbstractMatrix::operator*(const AbstractMatrix & other) const + { + const enumerator b = 32; + assert(Cols() == other.Rows()); + Matrix::type> ret(Rows(), other.Cols()); //check RVO + for (enumerator i0 = 0; i0 < Rows(); i0+=b) //loop rows + { + for (enumerator j0 = 0; j0 < other.Cols(); j0 += b) //loop columns + { + enumerator i0e = std::min(i0+b,Rows()); + enumerator j0e = std::min(j0+b,other.Cols()); + for (enumerator i = i0; i < i0e; ++i) + for(enumerator j = j0; j < j0e; ++j) + ret(i,j) = 0.0; + for (enumerator k0 = 0; k0 < Cols(); k0+=b) + { + enumerator k0e = std::min(k0+b,Cols()); + for (enumerator i = i0; i < i0e; ++i) + { + for (enumerator k = k0; k < k0e; ++k) + { + for(enumerator j = j0; j < j0e; ++j) + assign(ret(i, j),ret(i, j)+(*this)(i, k)*other(k, j)); + } + } + } + } + } + return ret; + } + */ + template template AbstractMatrix & @@ -1594,6 +1632,15 @@ namespace INMOST { // for A^T B assert(Rows() == B.Rows()); + /* + if( Rows() != Cols() ) + { + Matrix At = this->Transpose(); //m by n matrix + return (At*(*this)).Solve(At*B,print_fail); + } + Matrix AtB = B; //m by l matrix + Matrix AtA = (*this); //m by m matrix + */ Matrix At = this->Transpose(); //m by n matrix Matrix::type> AtB = At*B; //m by l matrix Matrix AtA = At*(*this); //m by m matrix @@ -1603,66 +1650,75 @@ namespace INMOST enumerator * order = new enumerator [m]; std::pair::type>,bool> ret = std::make_pair(Matrix::type>(m,l),true); - Var max, temp; - typename Promote::type tempb; + //Var temp; + INMOST_DATA_REAL_TYPE max,v; + //typeB tempb; for(enumerator i = 0; i < m; ++i) order[i] = i; for(enumerator i = 0; i < m; i++) { - enumerator maxk = i, maxq = i, temp2; - max = fabs(AtA(maxk,maxq)); + enumerator maxk = i, maxq = i;//, temp2; + max = fabs(get_value(AtA(maxk,maxq))); //Find best pivot - for(enumerator k = i; k < m; k++) // over rows + if( max < 1.0e-8 ) { - for(enumerator q = i; q < m; q++) // over columns + for(enumerator k = i; k < m; k++) // over rows { - if( fabs(AtA(k,q)) > max ) + for(enumerator q = i; q < m; q++) // over columns { - max = fabs(AtA(k,q)); - maxk = k; - maxq = q; + v = fabs(get_value(AtA(k,q))); + if( v > max ) + { + max = v; + maxk = k; + maxq = q; + } } } - } - //Exchange rows - if( maxk != i ) - { - for(enumerator q = 0; q < m; q++) // over columns of A - { - temp = AtA(maxk,q); - AtA(maxk,q) = AtA(i,q); - AtA(i,q) = temp; - } - //exchange rhs - for(enumerator q = 0; q < l; q++) // over columns of B - { - tempb = AtB(maxk,q); - AtB(maxk,q) = AtB(i,q); - AtB(i,q) = tempb; - } - } - //Exchange columns - if( maxq != i ) - { - for(enumerator k = 0; k < m; k++) //over rows + //Exchange rows + if( maxk != i ) { - temp = AtA(k,maxq); - AtA(k,maxq) = AtA(k,i); - AtA(k,i) = temp; + for(enumerator q = 0; q < m; q++) // over columns of A + { + std::swap(AtA(maxk,q),AtA(i,q)); + //temp = AtA(maxk,q); + //AtA(maxk,q) = AtA(i,q); + //AtA(i,q) = temp; + } + //exchange rhs + for(enumerator q = 0; q < l; q++) // over columns of B + { + std::swap(AtB(maxk,q),AtB(i,q)); + //tempb = AtB(maxk,q); + //AtB(maxk,q) = AtB(i,q); + //AtB(i,q) = tempb; + } } - //remember order in sol + //Exchange columns + if( maxq != i ) { - temp2 = order[maxq]; - order[maxq] = order[i]; - order[i] = temp2; + for(enumerator k = 0; k < m; k++) //over rows + { + std::swap(AtA(k,maxq),AtA(k,i)); + //temp = AtA(k,maxq); + //AtA(k,maxq) = AtA(k,i); + //AtA(k,i) = temp; + } + //remember order in sol + { + std::swap(order[maxq],order[i]); + //temp2 = order[maxq]; + //order[maxq] = order[i]; + //order[i] = temp2; + } } } //Check small entry - if( fabs(AtA(i,i)) < 1.0e-54 ) + if( fabs(get_value(AtA(i,i))) < 1.0e-54 ) { bool ok = true; for(enumerator k = 0; k < l; k++) // over columns of B { - if( fabs(AtB(i,k)/1.0e-54) > 1 ) + if( fabs(get_value(AtB(i,k))/1.0e-54) > 1 ) { ok = false; break; @@ -1815,8 +1871,8 @@ namespace INMOST /// @param coef Constant coefficient multiplying matrix. /// @param other Matrix to be multiplied. /// @return Matrix, each entry multiplied by a constant. -template -INMOST::Matrix::type> operator *(INMOST_DATA_REAL_TYPE coef, const INMOST::Matrix & other) +template +INMOST::Matrix::type> operator *(INMOST_DATA_REAL_TYPE coef, const INMOST::AbstractMatrix & other) {return other*coef;} #if defined(USE_AUTODIFF) /// Multiplication of matrix by a variable from left. @@ -1824,8 +1880,8 @@ INMOST::Matrix::type> oper /// @param coef Variable coefficient multiplying matrix. /// @param other Matrix to be multiplied. /// @return Matrix, each entry multiplied by a variable. -template -INMOST::Matrix::type> operator *(const INMOST::variable & coef, const INMOST::Matrix & other) +template +INMOST::Matrix::type> operator *(const INMOST::variable & coef, const INMOST::AbstractMatrix & other) {return other*coef;} /// Multiplication of matrix by a variable with first and /// second order derivatives from left. @@ -1833,8 +1889,8 @@ INMOST::Matrix::type> operator /// @param coef Variable coefficient multiplying matrix. /// @param other Matrix to be multiplied. /// @return Matrix, each entry multiplied by a variable. -template -INMOST::Matrix::type> operator *(const INMOST::hessian_variable & coef, const INMOST::Matrix & other) +template +INMOST::Matrix::type> operator *(const INMOST::hessian_variable & coef, const INMOST::AbstractMatrix & other) {return other*coef;} #endif diff --git a/Source/Headers/inmost_expression.h b/Source/Headers/inmost_expression.h index 471df6c257f3a6e4f6f8fdef15bfc5e97fd81481..90fa082733436b5b265a9b39e81e51e32ae7bf03 100644 --- a/Source/Headers/inmost_expression.h +++ b/Source/Headers/inmost_expression.h @@ -339,6 +339,11 @@ namespace INMOST std::cout << value << std::endl; entries.Print(); } + void swap(multivar_expression & b) + { + std::swap(value,b.value); + entries.Swap(b.entries); + } friend class multivar_expression_reference; }; @@ -351,6 +356,12 @@ namespace INMOST Sparse::Row entries; Sparse::HessianRow hessian_entries; public: + void swap(hessian_multivar_expression & b) + { + std::swap(value,b.value); + entries.Swap(b.entries); + hessian_entries.Swap(b.hessian_entries); + } /// Sets zero value and no first or second order variations. hessian_multivar_expression() :value(0) {} /// Sets value and no first or second order variations. @@ -1869,7 +1880,10 @@ template __INLINE INMOST_DATA_REAL_TY __INLINE void set_value(INMOST::multivar_expression_reference & Arg, const INMOST::multivar_expression & Val) {Arg.SetValue(Val.GetValue()); } __INLINE void set_value(INMOST::multivar_expression_reference & Arg, const INMOST::multivar_expression_reference & Val) {Arg.SetValue(Val.GetValue()); } __INLINE void assign(INMOST::var_expression & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;} - __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;} + __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;} + __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, const INMOST::var_expression & Val) {Arg = Val.GetValue();} + __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, const INMOST::multivar_expression & Val) {Arg = Val.GetValue();} + __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, const INMOST::multivar_expression_reference & Val) {Arg = Val.GetValue();} __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, const INMOST::var_expression & Val) {Arg = Val.GetValue();} __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, const INMOST::multivar_expression & Val) {Arg = Val.GetValue();} __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, const INMOST::multivar_expression_reference & Val) {Arg = Val.GetValue();} @@ -1943,6 +1957,9 @@ template __INLINE INMOST::function_expression< #else //USE_AUTODIFF +__INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = Val;} +__INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;} +__INLINE void assign(INMOST_REAL_INTEGER_TYPE & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = Val;} __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;} __INLINE void set_value(INMOST_DATA_REAL_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;} __INLINE INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE Arg) {return Arg;} diff --git a/Source/Headers/inmost_mesh.h b/Source/Headers/inmost_mesh.h index 985b1ae010ecd4da1fd1f3f349d075384222456d..0c23e43bd52cb0c3de21c4cda63507d81f199b33 100644 --- a/Source/Headers/inmost_mesh.h +++ b/Source/Headers/inmost_mesh.h @@ -3097,6 +3097,10 @@ namespace INMOST // types for BARYCENTER: EDGE | FACE | CELL // types for NORMAL: FACE | CELL (may precompute normal for cells in 2d case) // types for ORIENTATION: FACE + /// Marks face with the orientation direction by marker. + /// If marker is set then face is reversed. + /// Then all faces are oriented either inside or outside of the cell. + void FacesOrientation(ElementArray & faces, MarkerType rev); void PrepareGeometricData(GeomParam table); void RemoveGeometricData(GeomParam table); bool HaveGeometricData (GeometricData type, ElementType mask) const {return remember[type][ElementNum(mask)-1];} // requests to only one geometric and element type allowed diff --git a/Source/IO/mesh_ecl_file.cpp b/Source/IO/mesh_ecl_file.cpp index a4953f30ee247c0d228f3ac4f0778ec4f435bde3..df2c89e585d8a94b91e90f6e68e158807790f840 100644 --- a/Source/IO/mesh_ecl_file.cpp +++ b/Source/IO/mesh_ecl_file.cpp @@ -586,15 +586,15 @@ namespace INMOST { if (dir == "X-" || dir == "I-") return 0; - else if (dir == "X" || dir == "I" || dir == "X+" || dir == "I+") + else if (dir == "X" || dir == "I" || dir == "X+" || dir == "I+" || dir == "X " || dir == "I ") return 1; else if (dir == "Y-" || dir == "J-") return 2; - else if (dir == "Y" || dir == "J" || dir == "Y+" || dir == "J+") + else if (dir == "Y" || dir == "J" || dir == "Y+" || dir == "J+" || dir == "Y " || dir == "J ") return 3; else if (dir == "Z-" || dir == "K-") return 4; - else if (dir == "Z" || dir == "K" || dir == "Z+" || dir == "K+") + else if (dir == "Z" || dir == "K" || dir == "Z+" || dir == "K+" || dir == "Z " || dir == "K ") return 5; return -1; } @@ -2614,7 +2614,7 @@ namespace INMOST if (fault_cur.dir == -1) { if (verbosity > 0) - std::cout << "Unknown direction for fault " << fault_cur.name << " " << rec << std::endl; + std::cout << "Unknown direction for fault " << fault_cur.name << " " << rec << " (size:" < 0) std::cout << std::endl; + if (verbosity > 0) + { + if(project_perm ) + std::cout << "Permeability is projected onto axis" << std::endl; + else + std::cout << "Permeability is original" << std::endl; + } #if defined(USE_OMP) #pragma omp parallel #endif diff --git a/Source/Mesh/geometry.cpp b/Source/Mesh/geometry.cpp index dc0d2d6f1dabe3520fdc4063165c5f66942bcae9..a174bedb43567599e76d8f53b873bd05a2fd4d2e 100644 --- a/Source/Mesh/geometry.cpp +++ b/Source/Mesh/geometry.cpp @@ -27,6 +27,13 @@ namespace INMOST vecout[i] = vecin1[i] - vecin2[i]; } + __INLINE static void vec_diff(const Storage::real_array & vecin1,const Storage::real_array & vecin2, Storage::real * vecout, unsigned int size) + { + for(unsigned int i = 0; i < size; i++) + vecout[i] = vecin1[i] - vecin2[i]; + } + + __INLINE static void vec_cross_product(const Storage::real * vecin1, const Storage::real * vecin2, Storage::real * vecout) { Storage::real temp[3]; @@ -720,11 +727,118 @@ namespace INMOST } } + void Mesh::FacesOrientation(ElementArray & faces, MarkerType rev) + { + //can copy orientation-independent algorithm from + //incident_matrix.hpp: incident_matrix::compute_measure + //assume mdim is of size 3 at most + if( !faces.empty() ) + { + //real was = *ret/3.0; + Face cur = faces[0]; + Mesh * mesh = faces.GetMeshLink(); + //firstly, have to figure out orientation of each face + //mark all faces, so that we can perform adjacency retrival + MarkerType mrk = mesh->CreatePrivateMarker(); + //MarkerType rev = mesh->CreatePrivateMarker(); //reverse orientation + faces.SetPrivateMarker(mrk); //0-th face orientation is default + cur->RemPrivateMarker(mrk); + Node n1,n2; //to retrive edge + bool reverse = false; //reverse orientation in considered face + std::deque< orient_face > stack; //edge and first node and face for visiting + ElementArray edges = cur->getEdges(); + do + { + //figure out starting node order + if( edges[0]->getBeg() == edges[1]->getBeg() || + edges[0]->getBeg() == edges[1]->getEnd() ) + { + n1 = edges[0]->getEnd(); + n2 = edges[0]->getBeg(); + } + else + { + n1 = edges[0]->getBeg(); + n2 = edges[0]->getEnd(); + } + //schedule unvisited adjacent faces + for(unsigned j = 0; j < edges.size(); j++) + { + //schedule face adjacent to considered edge + ElementArray adjacent = edges[j]->getFaces(mrk); + assert(adjacent.size() <= 1); + if( !adjacent.empty() ) + { + adjacent.RemPrivateMarker(mrk); + stack.push_back(orient_face(edges[j],reverse ? n2 : n1,adjacent[0])); + } + //update edge nodes + n1 = n2; //current end is new begin + //find new end + if( n2 == edges[(j+1)%edges.size()]->getBeg() ) + n2 = edges[(j+1)%edges.size()]->getEnd(); + else + n2 = edges[(j+1)%edges.size()]->getBeg(); + } + if( stack.empty() ) break; + //get entry from stack + orient_face r = stack.front(); + //remove face from stack + stack.pop_front(); + //retrive edges for new face + edges = r.face->getEdges(); + reverse = false; + //figure out starting node order + if( edges[0]->getBeg() == edges[1]->getBeg() || + edges[0]->getBeg() == edges[1]->getEnd() ) + { + n1 = edges[0]->getEnd(); + n2 = edges[0]->getBeg(); + } + else + { + n1 = edges[0]->getBeg(); + n2 = edges[0]->getEnd(); + } + //find out common edge orientation + for(unsigned j = 0; j < edges.size(); j++) + { + if( edges[j] == r.bridge ) //found the edge + { + //reverse ordering on this face + if( r.first == n1 ) + { + if( isPrivate(rev) ) + r.face->SetPrivateMarker(rev); + else + r.face->SetMarker(rev); + reverse = true; + } + break; + } + //update edge nodes + n1 = n2; //current end is new begin + //find new end + if( n2 == edges[(j+1)%edges.size()]->getBeg() ) + n2 = edges[(j+1)%edges.size()]->getEnd(); + else + n2 = edges[(j+1)%edges.size()]->getBeg(); + } + } while(true); + faces.RemPrivateMarker(mrk); + mesh->ReleasePrivateMarker(mrk); + //faces.RemPrivateMarker(rev); + } + } + void Mesh::GetGeometricData(HandleType e, GeometricData type, Storage::real * ret) { assert(e != InvalidHandle()); assert(ret != NULL); - assert(type == MEASURE || type == CENTROID || type == BARYCENTER || type == NORMAL); + assert(type == MEASURE || + type == CENTROID || + type == BARYCENTER || + type == NORMAL); ElementType etype = GetHandleElementType(e); integer edim = Element::GetGeometricDimension(GetGeometricType(e)); integer mdim = GetDimensions(); @@ -746,8 +860,10 @@ namespace INMOST ElementArray nodes = Element(this,e)->getNodes(); if( nodes.size() > 1 ) { - Storage::real c[3]; - vec_diff(nodes[0]->Coords().data(),nodes[1]->Coords().data(),c,mdim); + real c[3] = {0,0,0}; + real_array v0 = nodes[0].Coords(); + real_array v1 = nodes[1].Coords(); + vec_diff(v0.data(),v1.data(),c,mdim); *ret = vec_len(c,mdim); } else *ret = 0; @@ -759,20 +875,26 @@ namespace INMOST ElementArray nodes = Element(this,e)->getNodes(); if( nodes.size() > 2 ) { - real x[3] = {0,0,0}; - Storage::real_array x0 = nodes[0].Coords(); - for(ElementArray::size_type i = 1; i < nodes.size()-1; i++) + *ret = 0; + real nt[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0},n0[3] = {0,0,0}, ss; + real_array v0 = nodes[0].Coords(); + real_array v1 = nodes[1].Coords(); + real_array v2 = nodes[2].Coords(); + vec_diff(v1,v0,l1,mdim); + vec_diff(v2,v0,l2,mdim); + vec_cross_product(l1,l2,n0); + for(int i = 1; i < (int)nodes.size()-1; i++) { - Storage::real_array v1 = nodes[i].Coords(); - Storage::real_array v2 = nodes[i+1].Coords(); - if( mdim == 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]); - } - x[2] += (v1[0]-x0[0])*(v2[1]-x0[1]) - (v1[1]-x0[1])*(v2[0]-x0[0]); + v1 = nodes[i].Coords(); + v2 = nodes[i+1].Coords(); + vec_diff(v1,v0,l1,mdim); + vec_diff(v2,v0,l2,mdim); + vec_cross_product(l1,l2,nt); + ss = vec_dot_product(n0,nt,3); + ss /= fabs(ss); + *ret += sqrt(vec_dot_product(nt,nt,3))*0.5*ss; } - *ret = ::sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])*0.5; + *ret = fabs(*ret); } else *ret = 0; //~ if( isnan(*ret) || fabs(*ret) < 1e-15 ) throw -1; break; @@ -782,133 +904,63 @@ namespace INMOST Cell me = Cell(this,e); ElementArray faces = me->getFaces(); - *ret = 0; - //can copy orientation-independent algorithm from - //incident_matrix.hpp: incident_matrix::compute_measure - //assume mdim is of size 3 at most - dynarray fcnt(mdim), fnrm(mdim), ccnt(mdim); - /* - me->Centroid(ccnt.data()); - for(ElementArray::size_type i = 0; i < faces.size(); i++) + bool ornt = !HaveGeometricData(ORIENTATION,FACE); + MarkerType rev = 0; + if( ornt ) { - faces[i]->Centroid(fcnt.data()); - for(int r = 0; r < mdim; ++r) - fcnt[r] = fcnt[r]-ccnt[r]; - faces[i]->OrientedNormal(me,fnrm.data()); - *ret += vec_dot_product(fcnt.data(),fnrm.data(),mdim); + rev = CreatePrivateMarker(); + FacesOrientation(faces,rev); } - if( *ret < 0.0 ) //a robust algorithm that can handle unoriented cell - */ - if( !faces.empty() ) + real vol = 0, a, at; + real x[3] = {0,0,0}, n[3] = {0,0,0}, n0[3] = {0,0,0}, s, ss; + real l1[3] = {0,0,0}, l2[3] = {0,0,0}; + real nt[3] = {0,0,0}; + for(unsigned j = 0; j < faces.size(); j++) { - //real was = *ret/3.0; - *ret = 0; - Face cur = faces[0]; - Cell c1 = me; - Mesh * mesh = c1.GetMeshLink(); - //firstly, have to figure out orientation of each face - //mark all faces, so that we can perform adjacency retrival - MarkerType mrk = mesh->CreatePrivateMarker(); - MarkerType rev = mesh->CreatePrivateMarker(); //reverse orientation - faces.SetPrivateMarker(mrk); //0-th face orientation is default - cur->RemPrivateMarker(mrk); - Node n1,n2; //to retrive edge - bool reverse = false; //reverse orientation in considered face - std::deque< orient_face > stack; //edge and first node and face for visiting - ElementArray edges = cur->getEdges(); - do + //compute normal to face + ElementArray nodes = faces[j].getNodes(); + if( ornt ) + s = faces[j].GetPrivateMarker(rev) ? -1.0 : 1.0; + else + s = faces[j].FaceOrientedOutside(me) ? 1.0 : -1.0; + x[0] = x[1] = x[2] = 0; + n[0] = n[1] = n[2] = 0; + a = 0; + real_array v0 = nodes[0].Coords(); + real_array v1 = nodes[1].Coords(); + real_array v2 = nodes[2].Coords(); + vec_diff(v1,v0,l1,mdim); + vec_diff(v2,v0,l2,mdim); + vec_cross_product(l1,l2,n0); + for(int k = 1; k < (int)nodes.size()-1; k++) { - //figure out starting node order - if( edges[0]->getBeg() == edges[1]->getBeg() || - edges[0]->getBeg() == edges[1]->getEnd() ) - { - n1 = edges[0]->getEnd(); - n2 = edges[0]->getBeg(); - } - else - { - n1 = edges[0]->getBeg(); - n2 = edges[0]->getEnd(); - } - //schedule unvisited adjacent faces - for(unsigned j = 0; j < edges.size(); j++) - { - //schedule face adjacent to considered edge - ElementArray adjacent = edges[j]->getFaces(mrk); - assert(adjacent.size() <= 1); - if( !adjacent.empty() ) - { - adjacent.RemPrivateMarker(mrk); - stack.push_back(orient_face(edges[j],reverse ? n2 : n1,adjacent[0])); - } - //update edge nodes - n1 = n2; //current end is new begin - //find new end - if( n2 == edges[(j+1)%edges.size()]->getBeg() ) - n2 = edges[(j+1)%edges.size()]->getEnd(); - else - n2 = edges[(j+1)%edges.size()]->getBeg(); - } - if( stack.empty() ) break; - //get entry from stack - orient_face r = stack.front(); - //remove face from stack - stack.pop_front(); - //retrive edges for new face - edges = r.face->getEdges(); - reverse = false; - //figure out starting node order - if( edges[0]->getBeg() == edges[1]->getBeg() || - edges[0]->getBeg() == edges[1]->getEnd() ) - { - n1 = edges[0]->getEnd(); - n2 = edges[0]->getBeg(); - } - else - { - n1 = edges[0]->getBeg(); - n2 = edges[0]->getEnd(); - } - //find out common edge orientation - for(unsigned j = 0; j < edges.size(); j++) - { - if( edges[j] == r.bridge ) //found the edge - { - //reverse ordering on this face - if( r.first == n1 ) - { - r.face->SetPrivateMarker(rev); - reverse = true; - } - break; - } - //update edge nodes - n1 = n2; //current end is new begin - //find new end - if( n2 == edges[(j+1)%edges.size()]->getBeg() ) - n2 = edges[(j+1)%edges.size()]->getEnd(); - else - n2 = edges[(j+1)%edges.size()]->getBeg(); - } - } while(true); - faces.RemPrivateMarker(mrk); - mesh->ReleasePrivateMarker(mrk); - c1->Centroid(ccnt.data()); - for(unsigned j = 0; j < faces.size(); j++) - { - //compute normal to face - faces[j].Centroid(fcnt.data()); - faces[j].Normal(fnrm.data()); - for(int r = 0; r < mdim; ++r) - fcnt[r] = fcnt[r]-ccnt[r]; - *ret += (faces[j]->GetPrivateMarker(rev) ? -1.0 : 1.0)*vec_dot_product(fcnt.data(),fnrm.data(),3); + v1 = nodes[k+0].Coords(); + v2 = nodes[k+1].Coords(); + vec_diff(v1,v0,l1,mdim); + vec_diff(v2,v0,l2,mdim); + vec_cross_product(l1,l2,nt); + for(int q = 0; q < 3; ++q) nt[q] *= 0.5; + ss = vec_dot_product(n0,nt,3); + ss /= fabs(ss); + at = sqrt(vec_dot_product(nt,nt,3))*ss; + //same as faces[j].Normal(n) + for(int q = 0; q < 3; ++q) + n[q] += nt[q]; + //same as faces[j].Centroid(x) + for(int q = 0; q < 3; ++q) + x[q] += at*(v0[q]+v1[q]+v2[q])/3.0; + a += at; } - *ret = fabs(*ret); + for(int q = 0; q < 3; ++q) x[q] /= a; + vol += s*vec_dot_product(x,n,3); + } + if( ornt ) + { + if( vol < 0.0 ) vol = -vol; faces.RemPrivateMarker(rev); - mesh->ReleasePrivateMarker(rev); + ReleasePrivateMarker(rev); } - *ret /= 3.0; - + *ret = vol/3.0; break; } } @@ -916,33 +968,32 @@ namespace INMOST //~ if( isnan(*ret) || fabs(*ret) < 1e-15 ) throw -1; break; case CENTROID: - if(etype == NODE ) - memcpy(ret,MGetDenseLink(e,CoordsTag()),sizeof(real)*mdim); - else if(HaveGeometricData(CENTROID,etype)) - { - memcpy(ret,MGetDenseLink(e,centroid_tag),sizeof(real)*mdim); - } - else - { - ElementArray nodes = Element(this,e)->getNodes(); - memset(ret,0,sizeof(real)*mdim); - if(nodes.size() != 0) + if(etype == NODE ) + memcpy(ret,MGetDenseLink(e,CoordsTag()),sizeof(real)*mdim); + else if(HaveGeometricData(CENTROID,etype)) + { + memcpy(ret,MGetDenseLink(e,centroid_tag),sizeof(real)*mdim); + } + else { - Storage::real div = 1.0/nodes.size(); - for(ElementArray::size_type i = 0; i < nodes.size(); i++) + ElementArray nodes = Element(this,e).getNodes(); + memset(ret,0,sizeof(real)*mdim); + for(unsigned k = 0; k < nodes.size(); ++k) { - Storage::real_array c =nodes[i].Coords(); - for(integer j = 0; j < mdim; j++) ret[j] += c[j]; + for(int q = 0; q < mdim; ++q) + ret[q] += nodes[k].Coords()[q]; } - for(integer j = 0; j < mdim; j++) ret[j] *= div; + for(int q = 0; q < mdim; ++q) + ret[q] /= (real)nodes.size(); } - } break; case BARYCENTER: - if( etype == NODE ) + if(etype == NODE ) memcpy(ret,MGetDenseLink(e,CoordsTag()),sizeof(real)*mdim); else if(HaveGeometricData(BARYCENTER,etype)) + { memcpy(ret,MGetDenseLink(e,barycenter_tag),sizeof(real)*mdim); + } else { memset(ret,0,sizeof(real)*mdim); @@ -951,112 +1002,132 @@ namespace INMOST ElementArray n = Element(this,e)->getNodes(); if( n.size() == 2 ) { - Storage::real_array a = n[0].Coords(); - Storage::real_array b = n[1].Coords(); - for(integer j = 0; j < dim; j++) - ret[j] = (a[j] + b[j])*0.5; + real_array v0 = n[0].Coords(); + real_array v1 = n[1].Coords(); + for(integer j = 0; j < dim; j++) + ret[j] = (v0[j] + v1[j])*0.5; } else if( n.size() == 1 ) { - Storage::real_array a = n[0].Coords(); - for(integer j = 0; j < dim; j++) ret[j] = a[j]; + real_array v0 = n[0].Coords(); + for(integer j = 0; j < dim; j++) ret[j] = v0[j]; } } else if( edim == 2 ) { ElementArray nodes = Element(this,e)->getNodes(); - real s,d, x1[3] = {0,0,0},x2[3] = {0,0,0},x[3] = {0,0,0}; - //here we compute area of polygon - //~ if( HaveGeometricData(MEASURE,etype) && HaveGeometricData(NORMAL,etype) ) - //~ { - //~ s = e->RealDF(measure_tag); - //~ memcpy(x,&e->RealDF(normal_tag),sizeof(Storage::real)*mdim); - //~ } - //~ else - //~ { if( nodes.size() > 2 ) { - Storage::real_array x0 = nodes[0].Coords(); - for(ElementArray::size_type 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 ) - { - 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]); - } - x[2] += (v1[0]-x0[0])*(v2[1]-x0[1]) - (v1[1]-x0[1])*(v2[0]-x0[0]); - } - s = ::sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]); - x[0] /= s; x[1] /= s; x[2] /= s; //here we obtain the unit normal - //~ } - //here we compute the center - Storage::real_array v0 = nodes[0].Coords(); - for(ElementArray::size_type j = 1; j < nodes.size()-1; j++) + *ret = 0; + real nt[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0}; + real c[3] = {0,0,0}, n0[3] = {0,0,0}, ss; + real_array v0 = nodes[0].Coords(); + real_array v1 = nodes[1].Coords(); + real_array v2 = nodes[2].Coords(); + vec_diff(v1,v0,l1,mdim); + vec_diff(v2,v0,l2,mdim); + vec_cross_product(l1,l2,n0); + real a = 0, at; + for(int i = 1; i < (int)nodes.size()-1; i++) { - Storage::real_array v1 = nodes[j].Coords(); - Storage::real_array v2 = nodes[j+1].Coords(); - for(integer k = 0; k < mdim; k++) - { - x1[k] = v0[k] - v1[k]; - x2[k] = v0[k] - v2[k]; - } - d = det3v(x1,x2,x); //here we use unit normal - for(integer k = 0; k < mdim; k++) - ret[k] += d*(v0[k]+v1[k]+v2[k]); + real_array v1 = nodes[i].Coords(); + real_array v2 = nodes[i+1].Coords(); + vec_diff(v1,v0,l1,mdim); + vec_diff(v2,v0,l2,mdim); + vec_cross_product(l1,l2,nt); + ss = vec_dot_product(n0,nt,3); + ss /= fabs(ss); + at = sqrt(vec_dot_product(nt,nt,3))*0.5*ss; + for(int q = 0; q < mdim; ++q) + c[q] += at*(v0[q]+v1[q]+v2[q])/3.0; + a += at; } - for(integer k = 0; k < mdim; k++) ret[k] /= 3.0 * s; + for(int q = 0; q < mdim; ++q) ret[q] = c[q]/a; } - else + //std::cout << ret[0] << " " << ret[1] << " " << ret[2] << std::endl; + } + else if( edim == 3 ) + { + Cell me = Cell(this,e); + ElementArray faces = me->getFaces(); + bool ornt = !HaveGeometricData(ORIENTATION,FACE); + MarkerType rev = 0; + if( ornt ) { - for(ElementArray::size_type i = 0; i < nodes.size(); i++) + rev = CreatePrivateMarker(); + FacesOrientation(faces,rev); + } + real vol = 0, a, at; + real x[3] = {0,0,0}, nt[3] = {0,0,0}, s; + real c[3] = {0,0,0}, n[3] = {0,0,0}; + real n0[3] = {0,0,0}, ss; + real l1[3] = {0,0,0}, l2[3] = {0,0,0}; + for(unsigned j = 0; j < faces.size(); j++) + { + //compute normal to face + ElementArray nodes = faces[j].getNodes(); + if( ornt ) + s = faces[j].GetPrivateMarker(rev) ? -1.0 : 1.0; + else + s = faces[j].FaceOrientedOutside(me) ? 1.0 : -1.0; + real_array v0 = nodes[0].Coords(); + real_array v1 = nodes[1].Coords(); + real_array v2 = nodes[2].Coords(); + vec_diff(v1,v0,l1,mdim); + vec_diff(v2,v0,l2,mdim); + vec_cross_product(l1,l2,n0); + x[0] = x[1] = x[2] = 0; + n[0] = n[1] = n[2] = 0; + a = 0; + for(int i = 1; i < (int)nodes.size()-1; i++) { - Storage::real_array c = nodes[i].Coords(); - for(integer k = 0; k < mdim; k++) ret[k] += c[k]; + real_array v1 = nodes[i].Coords(); + real_array v2 = nodes[i+1].Coords(); + vec_diff(v1,v0,l1,mdim); + vec_diff(v2,v0,l2,mdim); + vec_cross_product(l1,l2,nt); + for(int q = 0; q < 3; ++q) + nt[q] *= 0.5; + ss = vec_dot_product(n0,nt,3); + ss /= fabs(ss); + at = sqrt(vec_dot_product(nt,nt,3))*ss; + //same as faces[j].Normal(n) + for(int q = 0; q < 3; ++q) + n[q] += nt[q]; + //same as faces[j].Centroid(x) + for(int q = 0; q < 3; ++q) + x[q] += at*(v0[q]+v1[q]+v2[q])/3.0; + a += at; + //second-order midpoint formula + for(int q = 0; q < 3; ++q) + c[q] += s*nt[q]*(pow(v0[q]+v1[q],2)+pow(v0[q]+v2[q],2)+pow(v1[q]+v2[q],2))/24.0; } - for(integer k = 0; k < mdim; k++) ret[k] /= static_cast(nodes.size()); + for(int q = 0; q < 3; ++q) x[q] /= a; + vol += s*vec_dot_product(x,n,3); } - } - else if( edim == 3 ) - { - ElementArray faces = Element(this,e)->getFaces(); - real d,c,vol = 0, y[3]; - for(ElementArray::size_type i = 0; i < faces.size(); i++) + if( ornt ) { - d = y[0] = y[1] = y[2] = 0; - ElementArray nodes = faces[i].getNodes(); - Storage::real_array v0 = nodes[0].Coords(); - for(ElementArray::size_type j = 1; j < nodes.size()-1; j++) + if( vol < 0.0 ) { - Storage::real_array v1 = nodes[j].Coords(); - Storage::real_array v2 = nodes[j+1].Coords(); - c = det3v(v0.data(),v1.data(),v2.data()); - d += c; - y[0] += c * (v0[0] + v1[0] + v2[0]); - y[1] += c * (v0[1] + v1[1] + v2[1]); - y[2] += c * (v0[2] + v1[2] + v2[2]); + vol = -vol; + for(int q = 0; q < 3; ++q) + c[q] = -c[q]; } - c = faces[i].FaceOrientedOutside(Cell(this,e)) ? 1 : -1; - ret[0] += c * y[0]; - ret[1] += c * y[1]; - ret[2] += c * y[2]; - vol += c*d; + faces.RemPrivateMarker(rev); + ReleasePrivateMarker(rev); } - ret[0] /= vol*4; - ret[1] /= vol*4; - ret[2] /= vol*4; + vol /= 3.0; + for(int q = 0; q < mdim; ++q) + ret[q] = c[q]/(vol); + //std::cout << ret[0] << " " << ret[1] << " " << ret[2] << std::endl; } } break; case NORMAL: { -// real sret[3]; -// bool cmp = false; if( HaveGeometricData(NORMAL,etype) ) { memcpy(ret,MGetDenseLink(e,normal_tag),sizeof(real)*mdim); -// cmp = true; } else { @@ -1064,29 +1135,21 @@ namespace INMOST if( edim == 2 )//&& mdim == 3) { ElementArray nodes = Element(this,e)->getNodes(); - - Storage::real_array x0 = nodes[0].Coords(), a = x0, b; - for(ElementArray::size_type i = 0; i < nodes.size(); i++) + real n[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0}; + real nt[3] = {0,0,0}; + real_array v0 = nodes[0].Coords(); + for(int i = 1; i < (int)nodes.size()-1; i++) { - b = nodes[(i+1)%nodes.size()].Coords(); - ret[0] += (a[1]-x0[1])*(b[2]-x0[2]) - (a[2]-x0[2])*(b[1]-x0[1]); - ret[1] += (a[2]-x0[2])*(b[0]-x0[0]) - (a[0]-x0[0])*(b[2]-x0[2]); - ret[2] += (a[0]-x0[0])*(b[1]-x0[1]) - (a[1]-x0[1])*(b[0]-x0[0]); - a.swap(b); + real_array v1 = nodes[i].Coords(); + real_array v2 = nodes[i+1].Coords(); + vec_diff(v1,v0,l1,mdim); + vec_diff(v2,v0,l2,mdim); + vec_cross_product(l1,l2,nt); + for(int q = 0; q < 3; ++q) + n[q] += nt[q]*0.5; } - /* - for(unsigned i = 0; i < nodes.size(); i++) - { - Storage::real_array a = nodes[i].Coords(); - Storage::real_array b = nodes[(i+1)%nodes.size()].Coords(); - ret[0] += a[1]*b[2] - a[2]*b[1]; - ret[1] += a[2]*b[0] - a[0]*b[2]; - ret[2] += a[0]*b[1] - a[1]*b[0]; - } - */ - ret[0] *= 0.5; - ret[1] *= 0.5; - ret[2] *= 0.5; + for(int q = 0; q < mdim; ++q) + ret[q] = n[q]; } else if( edim == 1 )//&& mdim == 2 ) { @@ -1108,16 +1171,6 @@ namespace INMOST ret[1] *= l; } } -// real err = 0; -// for(int k = 0; k < 3; ++k) -// err += (ret[0]-sret[0])*(ret[0]-sret[0]); -// err = sqrt(err); -// if( err > 1.0e-8 && cmp ) -// { -// std::cout << "face " << GetHandleID(e) << " rank " << GetProcessorRank() << std::endl; -// std::cout << "cn " << ret[0] << " " << ret[1] << " " << ret[2] << std::endl; -// std::cout << "sn " << sret[0] << " " << sret[1] << " " << sret[2] << std::endl; -// } } } break; @@ -1298,15 +1351,12 @@ namespace INMOST data.RemPrivateMarker(mrk); mesh->ReleasePrivateMarker(mrk); Storage::real nrm[3], cnt[3], ccnt[3]; - c1->Centroid(ccnt); for(unsigned j = 0; j < data.size(); j++) { //std::cout << (data[j].GetPrivateMarker(rev) ? 0:1); //compute normal to face data[j].Centroid(cnt); data[j].Normal(nrm); - for(int r = 0; r < 3; ++r) - cnt[r] = cnt[r]-ccnt[r]; measure += (data[j]->GetPrivateMarker(rev) ? -1.0 : 1.0)*vec_dot_product(cnt,nrm,3); } //std::cout << "cur" << (cur->GetPrivateMarker(rev) ? 0:1) << " " << measure << " "; @@ -1526,7 +1576,7 @@ namespace INMOST { case MEASURE: return measure_tag; case CENTROID: return centroid_tag; - case BARYCENTER: return barycenter_tag; + case BARYCENTER: return barycenter_tag; case NORMAL: return normal_tag; } assert(false); diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index 0d14603b93f68eeb055d7dedd0b6c9ebe0bc1563..08899c8d0473ba7e53ae814f0762e11e5acbbb4d 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -1,3 +1,7 @@ +if(USE_MESH) +add_subdirectory(geom_test000) +endif(USE_MESH) + if(USE_AUTODIFF) add_subdirectory(autodiff_test000) add_subdirectory(autodiff_test001) @@ -18,4 +22,6 @@ add_subdirectory(pmesh_test001) endif() endif() + + add_subdirectory(xml_reader_test000) diff --git a/Tests/geom_test000/CMakeLists.txt b/Tests/geom_test000/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..bd18f992768d5e31c13dfef538047ceafa2c5e16 --- /dev/null +++ b/Tests/geom_test000/CMakeLists.txt @@ -0,0 +1,18 @@ +project(geom_test000) +set(SOURCE main.cpp) + +add_executable(geom_test000 ${SOURCE}) +target_link_libraries(geom_test000 inmost) + +if(USE_MPI) + message("linking geom_test000 with MPI") + target_link_libraries(geom_test000 ${MPI_LIBRARIES}) + if(MPI_LINK_FLAGS) + set_target_properties(geom_test000 PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") + endif() +endif(USE_MPI) + +add_test(NAME geom_test000_dual4 COMMAND $ ${CMAKE_CURRENT_SOURCE_DIR}/d4.pmf) +add_test(NAME geom_test000_cube4 COMMAND $ ${CMAKE_CURRENT_SOURCE_DIR}/c4.pmf) +add_test(NAME geom_test000_tetra4 COMMAND $ ${CMAKE_CURRENT_SOURCE_DIR}/t4.pmf) +add_test(NAME geom_test000_dual4_split COMMAND $ ${CMAKE_CURRENT_SOURCE_DIR}/d4sp.pmf) diff --git a/Tests/geom_test000/c4.pmf b/Tests/geom_test000/c4.pmf new file mode 100644 index 0000000000000000000000000000000000000000..c82c735ab564f1987e26957803476f50a43243ab Binary files /dev/null and b/Tests/geom_test000/c4.pmf differ diff --git a/Tests/geom_test000/d4.pmf b/Tests/geom_test000/d4.pmf new file mode 100644 index 0000000000000000000000000000000000000000..8975ea15bd52f1701a2a439dc56246fe32347d04 Binary files /dev/null and b/Tests/geom_test000/d4.pmf differ diff --git a/Tests/geom_test000/d4sp.pmf b/Tests/geom_test000/d4sp.pmf new file mode 100644 index 0000000000000000000000000000000000000000..c195ff0c42a7b0299923928361bc886b45d20f30 Binary files /dev/null and b/Tests/geom_test000/d4sp.pmf differ diff --git a/Tests/geom_test000/main.cpp b/Tests/geom_test000/main.cpp new file mode 100644 index 0000000000000000000000000000000000000000..7fdb6c1dd3b93db680f9e096a5c16a6f076d775d --- /dev/null +++ b/Tests/geom_test000/main.cpp @@ -0,0 +1,424 @@ +#include +#include + +#include "inmost.h" +using namespace INMOST; + +typedef Storage::real_array real_array; +typedef Storage::real real; + +#define ORIGINAL + +__INLINE static void vec_cross_product(const Storage::real * vecin1, const Storage::real * vecin2, Storage::real * vecout) +{ + Storage::real temp[3]; + temp[0] = vecin1[1]*vecin2[2] - vecin1[2]*vecin2[1]; + temp[1] = vecin1[2]*vecin2[0] - vecin1[0]*vecin2[2]; + temp[2] = vecin1[0]*vecin2[1] - vecin1[1]*vecin2[0]; + vecout[0] = temp[0]; + vecout[1] = temp[1]; + vecout[2] = temp[2]; +} + + +int main(int argc,char ** argv) +{ + int errors = 0; + Mesh::Initialize(&argc,&argv); + { +#if defined(ORIGINAL) + std::cout << "inmost algorithms" << std::endl; +#else + std::cout << "non-inmost algorithms" << std::endl; +#endif + Mesh m; + m.Load((argc>1)?argv[1]:"d4.pmf"); + rMatrix c(3,1), n(3,1), x(3,1), y(3,1), z(3,1), ct(3,1), nt(3,1); + rMatrix e1(3,1,0.0), e2(3,1,0.0), e3(3,1,0.0); + rMatrix x0(3,1),y0(3,1),z0(3,1); + rMatrix Et(3,3), Ef(3,3), d(3,1), W(9,3); + rMatrix nsum(3,1), ntsum(3,1); + e1(0,0) = 1.0; + e2(1,0) = 1.0; + e3(2,0) = 1.0; + Mesh::GeomParam t; + t[BARYCENTER] = CELL|FACE; + t[ORIENTATION] = FACE; + t[NORMAL] = FACE; + m.RemoveGeometricData(t); + //m.PrepareGeometricData(t); + for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) + it->FixNormalOrientation(); + int np = 0; + for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) + if( !it->Planarity() ) np++; + std::cout << "nonplanar: " << np << std::endl; + + /* + for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) + { + ElementArray faces = it->getFaces(); + rMatrix R(faces.size(),3), N(faces.size(),3), N2(3,faces.size()), R2(3,faces.size()); + for(unsigned k = 0; k < faces.size(); ++k) + { + faces[k].Barycenter(c.data()); + faces[k].OrientedNormal(it->self(),n.data()); + R(k,k+1,0,3) = c.Transpose(); + N(k,k+1,0,3) = n.Transpose(); + } + std::cout << "inv: " << std::endl; + (N.Transpose()*R).Invert().first.Print(); + N2 = (N.Transpose()*R).Invert().first*N.Transpose(); + R2 = (R.Transpose()*N).Invert().first*R.Transpose(); + nsum.Zero(); + for(unsigned k = 0; k < faces.size(); ++k) + nsum += N2(0,3,k,k+1); + std::cout << "Original N: " << std::endl; + N.Transpose().Print(); + std::cout << "Modified N: " << std::endl; + (N2*it->Volume()).Print(); + std::cout << "Original R: " << std::endl; + R.Transpose().Print(); + std::cout << "Modified R: " << std::endl; + (R2*it->Volume()).Print(); + nsum.Transpose().Print(); + (N2*R).Print(); + //((N.Transpose()*R).Invert().first*N.Transpose()*R).Print(); + } + return 0; + */ + /* + for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) + { + ElementArray faces = it->getFaces(); + rMatrix E(3,3),NN(3,3), W(9,9),K(3,3); + rMatrix R(faces.size(),3), N(faces.size(),3), R2(3,faces.size()),RK(faces.size(),3); + real v = 0, vk = 0; + + E.Zero(); N.Zero(); W.Zero(); + for(ElementArray::iterator jt = faces.begin(); jt != faces.end(); ++jt) + { + jt->Barycenter(c.data()); + jt->OrientedNormal(it->self(),n.data()); + E += c*n.Transpose()-(c.DotProduct(n)/3.0*rMatrix::Unit(3)); + NN += n*n.Transpose(); + for(int q = 0; q < 3; ++q) + { + W(q*4,0) += n(0,0)*n(0,0)/3.0; + W(q*4,1) += n(1,0)*n(0,0)/3.0; + W(q*4,2) += n(2,0)*n(0,0)/3.0; + W(q*4,3) += n(0,0)*n(1,0)/3.0; + W(q*4,4) += n(1,0)*n(1,0)/3.0; + W(q*4,5) += n(2,0)*n(1,0)/3.0; + W(q*4,6) += n(0,0)*n(2,0)/3.0; + W(q*4,7) += n(1,0)*n(2,0)/3.0; + W(q*4,8) += n(2,0)*n(2,0)/3.0; + } + } + std::cout << "Cell " << it->LocalID() << " vol " << it->Volume() << std::endl; + std::cout << "E:" << std::endl; + E.Print(); + std::cout << "N:" << std::endl; + N.Print(); + W += rMatrix::Unit(3).Kronecker(NN); + std::cout << "W" << std::endl; + W.Print(); + K = W.Solve(E.Repack(9,1),true).first.Repack(3,3); + std::cout << "K" << std::endl; + K.Print(); + E.Zero(); N.Zero(); + int k = 0; + for(ElementArray::iterator jt = faces.begin(); jt != faces.end(); ++jt) + { + jt->Barycenter(c.data()); + jt->OrientedNormal(it->self(),n.data()); + R(k,k+1,0,3) = c.Transpose(); + c -= K*n; + E += c*n.Transpose()-c.DotProduct(n)/3.0*rMatrix::Unit(3); + N(k,k+1,0,3) = n.Transpose(); + RK(k,k+1,0,3) = c.Transpose(); + k++; + } + v = 0; vk = 0; + R2 = (R.Transpose()*N).Invert().first*R.Transpose(); + iMatrix nf(1,faces.size()); + for(int k = 0; k < faces.size(); ++k) + { + //c = R(0,3,k,k+1); + c = R(k,k+1,0,3).Transpose(); + n = N(k,k+1,0,3).Transpose(); + v += c.DotProduct(n)/3.0; + + c = RK(k,k+1,0,3).Transpose(); + vk += c.DotProduct(n)/3.0; + + nf(0,k) = faces[k].LocalID(); + } + std::cout << "v " << v << " vk " << vk << std::endl; + std::cout << "E'" << std::endl; + E.Print(); + std::cout << "faces:" << std::endl; + nf.Print(); + std::cout << "R" << std::endl; + R.Transpose().Print(); + std::cout << "R2" << std::endl; + (R2*v).Print(); + std::cout << "||R2-R||: " << (R2*v-R).FrobeniusNorm() << std::endl; + std::cout << "RK" << std::endl; + RK.Transpose().Print(); + std::cout << "||RK-R||: " << (RK-R).FrobeniusNorm() << std::endl; + } + return 0; + */ + for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) + { + ElementArray faces = it->getFaces(); + //face center of mass, face normal, face area, cell volume + x.Zero(), y.Zero(), z.Zero(); + real vol = 0; + for(ElementArray::iterator jt = faces.begin(); jt != faces.end(); ++jt) + { +#if defined(ORIGINAL) + jt->Barycenter(c.data()); + jt->OrientedNormal(it->self(),n.data()); + x += c(0,0) * n; + y += c(1,0) * n; + z += c(2,0) * n; +#else + double s = jt->FaceOrientedOutside(it->self()) ? 1.0 : -1.0; + ElementArray nodes = jt->getNodes(); + real_array v0 = nodes[0].Coords(); + real_array v1 = nodes[1].Coords(); + real_array v2 = nodes[2].Coords(); + real a = 0, at; + n.Zero(); + c.Zero(); + real i1 = 0, i2 = 0, ss; + Et.Zero(); + rMatrix n0(3,1); + for(int q = 0; q < 3; ++q) + { + n0(q,0) = (v1[(q+1)%3]-v0[(q+1)%3])*(v2[(q+2)%3]-v0[(q+2)%3]); + n0(q,0)-= (v1[(q+2)%3]-v0[(q+2)%3])*(v2[(q+1)%3]-v0[(q+1)%3]); + n0(q,0)*=-1; + } + for(int k = 1; k < nodes.size()-1; ++k) + { + real_array v1 = nodes[k+0].Coords(); + real_array v2 = nodes[k+1].Coords(); + for(int q = 0; q < 3; ++q) + { + ct(q,0) = (v0[q]+v1[q]+v2[q])/3.0; + nt(q,0) = (v1[(q+1)%3]-v0[(q+1)%3])*(v2[(q+2)%3]-v0[(q+2)%3]); + nt(q,0)-= (v1[(q+2)%3]-v0[(q+2)%3])*(v2[(q+1)%3]-v0[(q+1)%3]); + } + ss = n0.DotProduct(nt); + ss /= fabs(ss); + nt *= 0.5; + i1 += ct.DotProduct(nt) * s; + //vol += ct.DotProduct(nt) * s; + at = sqrt(nt.DotProduct(nt)) * ss; + n += nt; + c += ct*at; + a += at; + Et += ct * nt.Transpose()*s; + //std::cout << "ct " << ct(0,0) << " " << ct(1,0) << " " << ct(2,0); + //std::cout << " nt " << nt(0,0) << " " << nt(1,0) << " " << nt(2,0); + //std::cout << " at " << at << std::endl; + //x += ct(0,0) * nt * s; + //y += ct(1,0) * nt * s; + //z += ct(2,0) * nt * s; + } + c /= a; + //std::cout << "c " << c(0,0) << " " << c(1,0) << " " << c(2,0); + //std::cout << " n " << n(0,0) << " " << n(1,0) << " " << n(2,0); + //std::cout << " a " << a << std::endl; + //i2 == i1! + Et -= i1/3.0*rMatrix::Unit(3); + + + Ef = c * n.Transpose()*s; + i2 = c.DotProduct(n) * s; + Ef -= i2/3.0*rMatrix::Unit(3); + + if( (Et-Ef).FrobeniusNorm() > 1.0e-9 ) + { + std::cout << "||Et-Ef||: " << (Et-Ef).FrobeniusNorm() << std::endl; + real_array v0 = nodes[0].Coords(); + real_array v1 = nodes[1].Coords(); + real_array v2 = nodes[2].Coords(); + n.Zero(); + a = 0; + rMatrix n0(3,1); + real l1[3],l2[3]; + for(int q = 0; q < 3; ++q) + { + l1[q] = v1[q]-v0[q]; + l2[q] = v2[q]-v0[q]; + } + vec_cross_product(l1,l2,n0.data()); + real ss; + for(int k = 1; k < nodes.size()-1; ++k) + { + v1 = nodes[k+0].Coords(); + v2 = nodes[k+1].Coords(); + for(int q = 0; q < 3; ++q) + { + ct(q,0) = (v0[q]+v1[q]+v2[q])/3.0; + //nt(q,0) = (v1[(q+1)%3]-v0[(q+1)%3])*(v2[(q+2)%3]-v0[(q+2)%3]); + //nt(q,0)-= (v1[(q+2)%3]-v0[(q+2)%3])*(v2[(q+1)%3]-v0[(q+1)%3]); + l1[q] = v1[q]-v0[q]; + l2[q] = v2[q]-v0[q]; + } + vec_cross_product(l1,l2,nt.data()); + ss = n0.DotProduct(nt); + ss /= fabs(ss); + nt *= 0.5; + n += nt; + at = sqrt(nt.DotProduct(nt))*ss; + a += at; + std::cout << "ct " << ct(0,0) << " " << ct(1,0) << " " << ct(2,0) << " "; + std::cout << "nt " << nt(0,0)/at << " " << nt(1,0)/at << " " << nt(2,0)/at << " "; + std::cout << "at " << at << std::endl; + } + real dd = n.DotProduct(rMatrix(v0.data(),3,1)); + std::cout << "n " << n(0,0)/a << " " << n(1,0)/a << " " << n(2,0)/a << " "; + std::cout << "planarity:"; + for(int k = 1; k < nodes.size(); ++k) + { + real_array v1 = nodes[k].Coords(); + std::cout << " " << n.DotProduct(rMatrix(v1.data(),3,1))-dd; + } + std::cout << std::endl; + // c = c+K*n + // Ef = c*n^T - c^T*n/3.0*I + // (Kc)*n^T - (Kc)^T*n/3.0*I = Et + // K*c*n^T - c^T*K^T*n/3.0*I = Et + // | k00 k01 k02 | nx + // (cx cy cz) | k10 k11 k12 | ny + // | k20 k21 k22 | nz + // d*n^T - d^T*n/3.0*I = Et + // + /* + W(0,0) = n(0,0) - n(0,0)/3.0; // e_xx + W(0,1) = - n(1,0)/3.0; // e_xx + W(0,2) = - n(2,0)/3.0; // e_xx + W(1,1) = n(1,0); // e_xy + W(2,2) = n(2,0); // e_xz + W(3,0) = n(0,0); // e_yx + W(4,0) = -n(0,0)/3.0; + W(4,1) = n(1,0) - n(1,0)/3.0; // e_yy + W(4,2) = -n(2,0)/3.0; + W(5,2) = n(2,0); // e_yz + W(6,0) = n(0,0); // e_zx + W(7,1) = n(1,0); // e_zy + W(8,0) = -n(0,0)/3.0; + W(8,1) = -n(1,0)/3.0; + W(8,2) = n(2,0)-n(2,0)/3.0; // e_zz + W *= s; + d = W.Solve((Et-Ef).Repack(9,1)).first; + c += d; + + Ef = c * n.Transpose()*s; + i2 = c.DotProduct(n) * s; + Ef -= i2/3.0*rMatrix::Unit(3); + + std::cout << "guess ||Et-Ef||: " << (Et-Ef).FrobeniusNorm() << std::endl; + */ + //std::cout << "Et: " << std::endl; + //Et.Print(); + //std::cout << "Before iters:" << std::endl; + //((c*n.Transpose() - (c.DotProduct(n)/3.0)*rMatrix::Unit(3))*s).Print(); + /* + vMatrix vn(3,1), vc(3,1), vE(9,1); + rMatrix W(9,6), E(9,1), up(6,1); + for(int nit = 0; nit < 50; ++nit) + { + for(int q = 0; q < 3; ++q) + { + vn(q,0) = unknown(n(q,0)+(rand()/(1.0*RAND_MAX))*1.0e-9,q); + vc(q,0) = unknown(c(q,0)+(rand()/(1.0*RAND_MAX))*1.0e-9,q); + } + vE = ((vc*vn.Transpose() - rMatrix::Unit(3)/3.0*vc.DotProduct(vn))*s - Et).Repack(9,1); + W.Zero(); + for(int q = 0; q < 9; ++q) + { + E(q,0) = vE(q,0).GetValue(); + for(Sparse::Row::iterator qt = vE(q,0).GetRow().Begin(); qt != vE(q,0).GetRow().End(); ++qt) + W(q,qt->first) = qt->second; + } + std::cout << nit << " " << E.DotProduct(E) << " "; + E.Transpose().Print(); + up = W.Solve(E).first; + n -= up(0,3,0,1); + c -= up(3,6,0,1); + } + + Ef = c * n.Transpose()*s; + i2 = c.DotProduct(n) * s; + Ef -= i2/3.0*rMatrix::Unit(3); + + std::cout << "iters ||Et-Ef||: " << (Et-Ef).FrobeniusNorm() << std::endl; + */ + + //std::cout << "After iters:" << std::endl; + //((c*n.Transpose() - (c.DotProduct(n)/3.0)*rMatrix::Unit(3))*s).Print(); + + } + + + + // + //std::cout << "area: " << a << std::endl; + + vol += c.DotProduct(n) * s; + x += c(0,0) * n * s; + y += c(1,0) * n * s; + z += c(2,0) * n * s; +#endif + } + + //std::cout << "nsum: " << nsum(0,0) << " " << nsum(1,0) << " " << nsum(2,0) << std::endl; + //std::cout << "ntsum: " << ntsum(0,0) << " " << ntsum(1,0) << " " << ntsum(2,0) << std::endl; + +#if defined(ORIGINAL) + x /= it->Volume(); + y /= it->Volume(); + z /= it->Volume(); + //std::cout << "inmost" << std::endl; + //x0.ConcatCols(y0).ConcatCols(z0).Print(); + //std::cout << "vol: " << it->Volume() << std::endl; +#else + vol /= 3.0; + x /= vol; + y /= vol; + z /= vol; + //std::cout << "non-inmost" << std::endl; + //x.ConcatCols(y).ConcatCols(z).Print(); + //std::cout << "vol: " << vol << std::endl; +#endif + if( (x-e1).FrobeniusNorm() > 1.0e-9 || x.CheckNans() ) + { + std::cout << "On CELL:" << it->LocalID() << " dx: " << x(0,0) << " " << x(1,0) << " " << x(2,0) << std::endl; + errors++; + } + if( (y-e2).FrobeniusNorm() > 1.0e-9|| y.CheckNans() ) + { + std::cout << "On CELL:" << it->LocalID() << " dy: " << y(0,0) << " " << y(1,0) << " " << y(2,0) << std::endl; + errors++; + } + if( (z-e3).FrobeniusNorm() > 1.0e-9|| z.CheckNans() ) + { + std::cout << "On CELL:" << it->LocalID() << " dz: " << z(0,0) << " " << z(1,0) << " " << z(2,0) << std::endl; + errors++; + } + } + } + Mesh::Finalize(); + if( errors ) + { + std::cout << "Errors: " << errors << std::endl; + return 1; + } + return 0; +} diff --git a/Tests/geom_test000/t4.pmf b/Tests/geom_test000/t4.pmf new file mode 100644 index 0000000000000000000000000000000000000000..711f3b6878fdfdd598e0b60b633ffcdc0825edab Binary files /dev/null and b/Tests/geom_test000/t4.pmf differ