MeshKit  1.0
SubMapping.cpp
Go to the documentation of this file.
00001 #include "meshkit/SubMapping.hpp"
00002 #include "meshkit/MKCore.hpp"
00003 #include "meshkit/iMesh.hpp"
00004 #include "meshkit/RegisterMeshOp.hpp"
00005 #include "meshkit/EdgeMesher.hpp"
00006 #include "meshkit/VertexMesher.hpp"
00007 #include "meshkit/TFIMapping.hpp"
00008 #include "meshkit/SimpleArray.hpp"
00009 #include "lp_lib.h"
00010 #include "meshkit/MeshImprove.hpp"
00011 #include "LPSolveClass.hpp"
00012 #include "IsoLaplace.hpp"
00013 #include "EquipotentialSmooth.hpp"
00014 #ifdef HAVE_MESQUITE
00015 #include "meshkit/MeshImprove.hpp"
00016 #endif
00017 #include <iostream>
00018 #include <algorithm>
00019 #include <math.h>
00020 #include <map>
00021 
00022 
00023 
00024 namespace MeshKit
00025 {
00026 
00027 //---------------------------------------------------------------------------//
00028 //Entity Type initilization for SubMapping meshing
00029 moab::EntityType SubMapping_tps[] = {moab::MBVERTEX, moab::MBEDGE, moab::MBQUAD, moab::MBMAXTYPE};
00030 const moab::EntityType* SubMapping::output_types()
00031   { return SubMapping_tps; }
00032 
00033 //---------------------------------------------------------------------------//
00034 // construction function for SubMapping class
00035 SubMapping::SubMapping(MKCore *mk_core, const MEntVector &me_vec) : MeshScheme(mk_core, me_vec)
00036 {
00037         //buildAssociation();
00038 }
00039 
00040 
00041 //---------------------------------------------------------------------------//
00042 // deconstruction function for SubMapping class
00043 SubMapping::~SubMapping()
00044 {
00045 
00046 }
00047 
00048 //---------------------------------------------------------------------------//
00049 // setup function: define the size between the different layers
00050 void SubMapping::setup_this()
00051 {
00052         if (mentSelection.empty())
00053         return;
00054     //compute the number of intervals for the associated ModelEnts, from the size set on them
00055     //the sizing function they point to, or a default sizing function
00056   for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
00057   {
00058       ModelEnt *me = mit -> first;
00059           int dimension = me->dimension();
00060           if (dimension != 2)
00061                  ECERRCHK(MK_FAILURE, "bad input for TFI Mapping, we can only mesh surfaces");
00062           //first check whether the surface is meshed or not
00063       if (me->get_meshed_state() >= COMPLETE_MESH)
00064            continue;
00065 
00066           
00067           SizingFunction *sf = mk_core()->sizing_function(me->sizing_function_index());
00068           if (!sf && me -> mesh_intervals() < 0 && me -> interval_firmness() == DEFAULT &&
00069         mk_core()->sizing_function(0))
00070         sf = mk_core()->sizing_function(0);
00071 
00072           if (!sf && me -> mesh_intervals() < 0 && me -> interval_firmness() == DEFAULT){
00073          //no sizing set, just assume default #intervals as 20
00074          me->mesh_intervals(20);
00075          me->interval_firmness(DEFAULT);
00076       }
00077       else{
00078                         //check # intervals first, then size, and just choose for now
00079                 if (sf->intervals() > 0)
00080                         throw Error(MK_INCOMPLETE_MESH_SPECIFICATION,  "Sizing function for edge should have edge size instead of intervals for submapping.");
00081                         else if (sf->size()>0){
00082                                 size_low_bound = sf->size();
00083                         me->interval_firmness(SOFT);
00084                         }
00085                         else
00086                         throw Error(MK_INCOMPLETE_MESH_SPECIFICATION,  "Sizing function for edge had neither positive size nor positive intervals.");
00087       }   
00088   }
00089 
00090 
00091 }
00092 
00093 
00094 //---------------------------------------------------------------------------//
00095 // execute function: generate the all-hex mesh through sweeping from source 
00096 // surface to target surface
00097 void SubMapping::execute_this()
00098 {
00099         std::vector<double> coords;
00100         std::vector<moab::EntityHandle> nodes;
00101 
00102         iBase_TagHandle global_tag, global_tag1;
00103         iGeom::Error g_err = mk_core()->igeom_instance()->getTagHandle("GLOBAL_ID", global_tag);
00104         IBERRCHK(g_err, "Trouble get the mesh entity set from geometric edges.");
00105 
00106         iMesh::Error m_err = mk_core()->igeom_instance()->getTagHandle("GLOBAL_ID", global_tag1);
00107         IBERRCHK(m_err, "Trouble get the mesh entity set from geometric edges.");
00108         
00109         for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
00110         {
00111         ModelEnt *me = mit -> first;            
00112                 //first check whether the surface is meshed or not
00113                 if (me->get_meshed_state() >= COMPLETE_MESH)
00114                         continue;
00115 
00116                 //classify the vertex as END, CORNER, REVERSAL, SIDE
00117                 VertexClassification(me);
00118         //classify the edge as +I, -I, +J and -J
00119                 EdgeClassification();
00120       //low bound for mesh size
00121                 EdgeDiscretization(me);
00122 
00123                 InteriorNodeInterpolation(me);
00124                 MeshSmoothing(me);
00125 
00126                 //ok, we are done!
00127                 vertices_types.clear();
00128                 nodes.clear();
00129                 vertices.clear();
00130                 edges.clear();
00131                 interior_angle.clear();
00132                 sorted_vertex_list.clear();
00133                 sorted_node_list.clear();
00134                 sorted_edge_list.clear();
00135                 edges_types.clear();
00136                 coordinate_i_j.clear();
00137                 edge_size.clear();
00138                 quads.clear();
00139 
00140                 me->commit_mesh(mit->second, COMPLETE_MESH);
00141                 mk_core()->save_mesh("test_cai.vtk");
00142         }
00143 
00144 }
00145 
00146 //setup the mesh size
00147 void SubMapping::SetupMeshSize(double size)
00148 {
00149         size_low_bound = size;
00150 }
00151 
00152 
00153 //classify the vertices as side, end, reversal and corner
00154 //end           --an interior angle of approximately 90 degrees
00155 //side          --an interior angle of approximately 180 degrees = 0
00156 //corner        --an interior angle of approximately 270 degrees = -90
00157 //reversal  --an interior angle of approximately 360 degrees = -180
00158 void SubMapping::VertexClassification(ModelEnt *ent)
00159 {
00160         //create a taghandle
00161         iGeom::Error g_err = mk_core()->igeom_instance()->getTagHandle("GEOM_SUBMAPPING", g_taghandle);
00162         if (g_err){
00163                 g_err = mk_core()->igeom_instance()->createTag("GEOM_SUBMAPPING", 1, iBase_INTEGER, g_taghandle);       
00164                 IBERRCHK(g_err, "Trouble create a taghandle.");
00165         }
00166 
00167         //get the geometric vertices from the surface
00168         vertices.clear();
00169         vector<iBase_EntityHandle> tmp;
00170         g_err = mk_core()->igeom_instance()->getEntAdj(ent->geom_handle(), iBase_VERTEX, tmp);
00171         IBERRCHK(g_err, "Trouble get the adjacent geometric vertices on a surface.");
00172         assert(tmp.size()>0);
00173 
00174         vertices.resize(tmp.size());
00175         for (unsigned int i = 0; i < tmp.size(); i++){
00176                 vertices[i].index = i;
00177                 vertices[i].gVertexHandle = tmp[i];
00178 
00179                 g_err = mk_core()->igeom_instance()->getVtxCoord(vertices[i].gVertexHandle, vertices[i].xyz[0], vertices[i].xyz[1], vertices[i].xyz[2]);
00180                 IBERRCHK(g_err, "Trouble get the coordinates from vertices.");
00181 
00182                 g_err = mk_core()->igeom_instance()->setIntData(vertices[i].gVertexHandle, g_taghandle, i);
00183                 IBERRCHK(g_err, "Trouble set the int data for the geometric vertex.");
00184         }
00185 
00186         vertices_types.resize(vertices.size());
00187         
00188         //get the geometric edges
00189         edges.clear();
00190         
00191         tmp.clear();
00192         g_err = mk_core()->igeom_instance()->getEntAdj(ent->geom_handle(), iBase_EDGE, tmp);
00193         IBERRCHK(g_err, "Trouble get the adjacent geometric vertices on a surface.");
00194 
00195         std::set<iBase_EntityHandle> edge_set;
00196         edges.resize(tmp.size());
00197         for (unsigned int i = 0; i < tmp.size(); i++){
00198                 edges[i].connect.resize(2);
00199                 edges[i].index = i;
00200                 edges[i].gEdgeHandle = tmp[i];
00201                 edge_set.insert(tmp[i]);
00202 
00203                 g_err = mk_core()->igeom_instance()->setIntData(edges[i].gEdgeHandle, g_taghandle, i);
00204                 IBERRCHK(g_err, "Trouble set the int data for the geometric edge.");
00205 
00206                 vector<iBase_EntityHandle> adj_vertices;
00207                 g_err = mk_core()->igeom_instance()->getEntAdj(edges[i].gEdgeHandle, iBase_VERTEX, adj_vertices);
00208                 IBERRCHK(g_err, "Trouble get the adjacent vertices for the geometric edge.");
00209 
00210                 assert(adj_vertices.size()<=2);
00211                 for (unsigned int j = 0; j < adj_vertices.size(); j++){
00212                         int index = -1;
00213                         g_err = mk_core()->igeom_instance()->getIntData(adj_vertices[j], g_taghandle, index);
00214                         IBERRCHK(g_err, "Trouble get the int data of geometric vertex.");
00215 
00216                         edges[i].connect[j] = &vertices[index];
00217                 }
00218         }
00219 
00220         //organize the vertices and edges on the boundaries
00221         VerEdgOrganize(edge_set, tmp, ent->geom_handle());
00222 
00223         //calculate the angle for vertices on the boundaries
00224         interior_angle.resize(vertices.size());
00225         GetAngle(ent->geom_handle(), interior_angle);
00226 
00227         //initial vertex classification,
00228         //need to use the linear programming method to produce the valid vertex classification 
00229         //if the initial vertex classification is not a valid vertex classification
00230         vertices_types.resize(vertices.size());
00231         for (unsigned int i = 0; i < interior_angle.size(); i++){
00232                 if ((interior_angle[i] >= 45)&&(interior_angle[i] <= 130)){
00233                         vertices_types[i] = END;
00234                 }
00235                 else if ((interior_angle[i] < 225)&&(interior_angle[i] > 130)){
00236                         vertices_types[i] = SIDE;
00237                 }
00238                 else if ((interior_angle[i] > 225)&&(interior_angle[i] < 315)){
00239                         vertices_types[i] = CORNER;
00240                 }
00241                 else{
00242                         vertices_types[i] = REVERSAL;
00243                 }
00244         }
00245 
00246         //extra step
00247         for (unsigned int m = 0; m < vertices_types.size(); m++){
00248                         
00249                 if (vertices_types[m] == REVERSAL)
00250                         vertices_types[m] = SIDE;
00251         }
00252 
00253 
00254 
00255 
00256         int sum = 0;
00257         for (std::vector<VertexTypes>::iterator it = vertices_types.begin(); it != vertices_types.end(); it++)
00258                 sum += *it;
00259 
00260         std::cout << "\n\n\nbefore, sum of vertex type is " << sum << std::endl;
00261         
00262         if (sum != 4){//check whether the initial vertex classification is valid or not
00263                 //use the linear programming to produce a valid vertex classification
00264                                 
00265                 VtxClassificationLP();
00266         }
00267         
00268 
00269         std::cout << "\n\n\nsum = " << sum << std::endl;
00270 }
00271 
00272 //use the lpsolve library to solve the linear programming
00273 void SubMapping::VtxClassificationLP()
00274 {
00275         LPSolveClass lp;
00276         vector<int> VtxType1(vertices_types.size()), VtxType2;
00277         int m = 0;
00278         for (std::vector<VertexTypes>::iterator it = vertices_types.begin(); it != vertices_types.end(); it++){
00279                         
00280                 if (*it == REVERSAL)
00281                         vertices_types[m] = SIDE;
00282                 
00283                 VtxType1[m] = int(vertices_types[m]);           
00284                 m++;
00285         }
00286 
00287         //setup the objective function: minimize the max value, first half variables are new variables, the latter half variables are old variables
00288         vector<double> b(2*vertices_types.size());
00289         for (unsigned int i = 0; i < vertices_types.size(); i++){
00290                 b[i] = 1.0;
00291                 b[vertices_types.size()+i] = 0.0;
00292         }
00293         lp.SetupObj(b, 0.0);
00294         
00295         //setup the equality constraint
00296         vector<vector<double> > left;
00297         vector<double> right;
00298         right.push_back(4.0);
00299         left.resize(1);
00300         left[0].resize(2*vertices_types.size());
00301         for (unsigned int i = 0; i < vertices_types.size(); i++){
00302                 left[0][i] = 0.0;
00303                 left[0][vertices_types.size()+i] = 1.0;
00304         }
00305         /*
00306         for (unsigned int i = 0; i < vertices_types.size(); i++){
00307                 right.push_back(0.0);
00308                 left[i+1].resize(2*vertices_types.size());
00309                 for (unsigned int j = 0; j < vertices_types.size(); j++){
00310                         left[i+1][j] = 0.0;
00311                         left[i+1][vertices_types.size()+j] = 0.0;
00312                         if ((i==j)&&(vertices_types[i]==END)){//fix the END vertex type
00313                                 left[i+1][vertices_types.size()+j] = 1.0;
00314                                 right[i+1] = 1.0;
00315                         }
00316                 }
00317         }
00318         */
00319         lp.SetupEqu(left, right);
00320 
00321         //setup the inequality constraint
00322         left.clear();
00323         right.clear();
00324         left.resize(6*vertices_types.size());
00325         right.resize(6*vertices_types.size());
00326         for (unsigned int i = 0; i < vertices_types.size(); i++){
00327                 right[i] = VtxType1[i];
00328                 right[vertices_types.size()+i] = VtxType1[i];
00329                 right[2*vertices_types.size()+i] = 1;
00330                 right[3*vertices_types.size()+i] = 2;
00331                 right[4*vertices_types.size()+i] = 1+VtxType1[i];
00332                 right[5*vertices_types.size()+i] = 1-VtxType1[i];
00333         }
00334         for (unsigned int i = 0; i < vertices_types.size(); i++){
00335                 left[i].resize(2*vertices_types.size());
00336                 left[vertices_types.size()+i].resize(2*vertices_types.size());
00337                 for (unsigned int j= 0; j < vertices_types.size(); j++){
00338                         left[i][j] = -1.0;
00339                         left[i][vertices_types.size()+j] = 1.0;
00340                         left[vertices_types.size()+i][j] = -1.0;
00341                         left[vertices_types.size()+i][vertices_types.size()+j] = -1.0;
00342                 }
00343                 //constrain the variable range
00344                 left[2*vertices_types.size()+i].resize(2*vertices_types.size());
00345                 left[3*vertices_types.size()+i].resize(2*vertices_types.size());
00346                 left[4*vertices_types.size()+i].resize(2*vertices_types.size());
00347                 left[5*vertices_types.size()+i].resize(2*vertices_types.size());
00348                 for (unsigned int j = 0; j < 2*vertices_types.size(); j++){
00349                         left[2*vertices_types.size()+i][j] = 0.0;
00350                         left[3*vertices_types.size()+i][j] = 0.0;
00351                         left[4*vertices_types.size()+i][j] = 0.0;
00352                         left[5*vertices_types.size()+i][j] = 0.0;
00353                 }
00354                 left[2*vertices_types.size()+i][vertices_types.size()+i] = 1.0;
00355                 left[3*vertices_types.size()+i][vertices_types.size()+i] = -1.0;
00356 
00357                 //constrain the difference between new vertex type and initial vertex type
00358                 left[4*vertices_types.size()+i][vertices_types.size()+i] = 1.0;
00359                 left[5*vertices_types.size()+i][vertices_types.size()+i] = -1.0;
00360                         
00361         }
00362         lp.SetupInEqu(left, right);
00363 
00364         lp.Execute();
00365         //get the solved variables
00366         vector<int> vtx_types;
00367         lp.GetVariables(vtx_types);
00368         for (unsigned int i = 0; i < vertices_types.size(); i++){
00369                 switch(vtx_types[vertices_types.size()+i]){
00370                         case 1:
00371                                 vertices_types[i] = END;
00372                                 break;
00373                         case 0:
00374                                 vertices_types[i] = SIDE;
00375                                 break;
00376                         case -1:
00377                                 vertices_types[i] = CORNER;
00378                                 break;
00379                         case -2:
00380                                 vertices_types[i] = REVERSAL;
00381                                 break;
00382                         default:
00383                                 break;
00384                 }
00385         }
00386 
00387         //check whether sum = 4;
00388         int sum = 0;
00389         for (std::vector<VertexTypes>::iterator it = vertices_types.begin(); it != vertices_types.end(); it++)
00390                 sum += *it;
00391         assert(sum == 4);
00392 }
00393 
00394 //calculate the angle for vertices on the boundaries
00395 void SubMapping::GetAngle(iBase_EntityHandle surf, vector<double> &angle)
00396 {
00397         angle.resize(sorted_vertex_list.size());
00398         for (unsigned int i = 0; i < sorted_vertex_list.size(); i++){
00399                 //vertices = i,  previous edge = (i+sorted_edge_list.size()-1)%sorted_edge_list.size(), next edge = (i+sorted_edge_list.size()+1)%sorted_edge_list.size()
00400 
00401                 //calculate the tangent vector at the previous edge
00402                 double tangvector_p[3];
00403                 int previous = (i+sorted_edge_list.size()-1)%sorted_edge_list.size();
00404                 //calculate the tangent vector and edge sense with respect to two vertices
00405                 iGeom::Error g_err = mk_core()->igeom_instance()->getEntTgntXYZ(edges[sorted_edge_list[previous]].gEdgeHandle, vertices[sorted_vertex_list[i]].xyz[0], vertices[sorted_vertex_list[i]].xyz[1], vertices[sorted_vertex_list[i]].xyz[2], tangvector_p[0], tangvector_p[1], tangvector_p[2]);
00406                 IBERRCHK(g_err, "Trouble get the tangent vector at a vertex of an edge.");
00407 
00408                 int sense = -10;
00409                 g_err = mk_core()->igeom_instance()->getEgVtxSense(edges[sorted_edge_list[previous]].gEdgeHandle, vertices[sorted_vertex_list[i]].gVertexHandle, vertices[sorted_vertex_list[(i+vertices.size()-1)%vertices.size()]].gVertexHandle, sense);
00410                 IBERRCHK(g_err, "Trouble get edge sense with respect to two vertices.");
00411 
00412                 int pre_edge_sense_face = -10;
00413                 g_err = mk_core()->igeom_instance()->getEgFcSense(edges[sorted_edge_list[previous]].gEdgeHandle, surf, pre_edge_sense_face);
00414                 IBERRCHK(g_err, "Trouble get edge sense with respect to the face.");
00415                 std::cout << "-------------------\nvertex index = " << i << "\tprevious edge sense with respect to face = " << pre_edge_sense_face << std::endl;
00416 
00417                 if (sense < 0){
00418                         tangvector_p[0] = -1.0*tangvector_p[0];
00419                         tangvector_p[1] = -1.0*tangvector_p[1];
00420                         tangvector_p[2] = -1.0*tangvector_p[2];
00421                 }
00422                                         
00423                 //calculate the tangent vector at next edge
00424                 double tangvector_n[3];
00425                 int next = (i+sorted_edge_list.size())%sorted_edge_list.size();
00426                 g_err = mk_core()->igeom_instance()->getEntTgntXYZ(edges[sorted_edge_list[next]].gEdgeHandle, vertices[sorted_vertex_list[i]].xyz[0], vertices[sorted_vertex_list[i]].xyz[1], vertices[sorted_vertex_list[i]].xyz[2], tangvector_n[0], tangvector_n[1], tangvector_n[2]);
00427                 IBERRCHK(g_err, "Trouble get the tangent vector at a vertex of an edge.");
00428                 g_err = mk_core()->igeom_instance()->getEgVtxSense(edges[sorted_edge_list[next]].gEdgeHandle, vertices[sorted_vertex_list[i]].gVertexHandle, vertices[sorted_vertex_list[(i+1)%vertices.size()]].gVertexHandle, sense);
00429                 IBERRCHK(g_err, "Trouble get edge sense with respect to two vertices.");
00430                 
00431                 int next_edge_sense_face = -10;
00432                 g_err = mk_core()->igeom_instance()->getEgFcSense(edges[sorted_edge_list[next]].gEdgeHandle, surf, next_edge_sense_face);
00433                 IBERRCHK(g_err, "Trouble get edge sense with respect to the face.");
00434                 std::cout << "\t" << i << "\tnext edge sense with respect to face = " << next_edge_sense_face << std::endl;
00435 
00436                 if (sense < 0){
00437                         tangvector_n[0] = -1.0*tangvector_n[0];
00438                         tangvector_n[1] = -1.0*tangvector_n[1];
00439                         tangvector_n[2] = -1.0*tangvector_n[2];
00440                 }
00441 
00442                 
00443                 //calculate the crossproduct
00444                 double crossproduct[3];
00445                 crossproduct[0] = tangvector_n[1]*tangvector_p[2]-tangvector_n[2]*tangvector_p[1];
00446                 crossproduct[1] = tangvector_n[2]*tangvector_p[0]-tangvector_n[0]*tangvector_p[2];
00447                 crossproduct[2] = tangvector_n[0]*tangvector_p[1]-tangvector_n[1]*tangvector_p[0];
00448 
00449                 //calculate the dotproduct
00450                 double dotproduct = tangvector_n[0]*tangvector_p[0] + tangvector_n[1]*tangvector_p[1] + tangvector_n[2]*tangvector_p[2];
00451                 double theta_pi = acos(dotproduct/(sqrt(pow(tangvector_n[0],2)+pow(tangvector_n[1], 2)+pow(tangvector_n[2], 2))*sqrt(pow(tangvector_p[0],2)+pow(tangvector_p[1],2)+pow(tangvector_p[2],2))));
00452                 double theta = theta_pi*180.0/3.1415926;
00453 
00454                 //get the surface norm at a specific point
00455                 double normal[3];
00456                 g_err = mk_core()->igeom_instance()->getEntNrmlXYZ(surf, vertices[sorted_vertex_list[i]].xyz[0], vertices[sorted_vertex_list[i]].xyz[1], vertices[sorted_vertex_list[i]].xyz[2], normal[0], normal[1], normal[2]);
00457                 IBERRCHK(g_err, "Trouble get the normal at a specific point on the surface.");
00458                 dotproduct = crossproduct[0]*normal[0]+crossproduct[1]*normal[1]+crossproduct[2]*normal[2];
00459                 
00460                 if (dotproduct < 0)
00461                         theta = 360.0 - theta;
00462                 
00463                 angle[sorted_vertex_list[i]] = theta;
00464         }
00465 
00466 }
00467 
00468 bool SubMapping::isCurved(int vtx_index, vector<double> u1, vector<double> u2, vector<double> u3, vector<double> u4, vector<vector<double> > tang_pre, vector<vector<double> > tang_next)
00469 {
00470         int previous = (vtx_index+sorted_edge_list.size()-1)%sorted_edge_list.size(), next = (vtx_index+sorted_edge_list.size())%sorted_edge_list.size();
00471         double  dotproduct, vec[3], theta1, theta2, pts2[3], pts1[3], u, length;
00472         u = u1[vtx_index] + 0.5*(u2[vtx_index]-u1[vtx_index]);
00473         iGeom::Error g_err = mk_core()->igeom_instance()->getEntUtoXYZ(edges[sorted_edge_list[previous]].gEdgeHandle, u, pts1[0], pts1[1], pts1[2]);
00474         IBERRCHK(g_err, "Trouble get xyz coordinates from parametric coordinates on edge.");
00475         u = u1[vtx_index] + 0.501*(u2[vtx_index]-u1[vtx_index]);
00476         g_err = mk_core()->igeom_instance()->getEntUtoXYZ(edges[sorted_edge_list[previous]].gEdgeHandle, u, pts2[0], pts2[1], pts2[2]);
00477         IBERRCHK(g_err, "Trouble get xyz coordinates from parametric coordinates on edge.");
00478         for (int j = 0; j < 3; j++)
00479                 vec[j] = pts2[j] - pts1[j];
00480         dotproduct = vec[0]*tang_pre[vtx_index][0] + vec[1]*tang_pre[vtx_index][1] + vec[2]*tang_pre[vtx_index][2];
00481         length = sqrt(pow(vec[0],2)+pow(vec[1], 2)+pow(vec[2], 2))*sqrt(pow(tang_pre[vtx_index][0],2)+pow(tang_pre[vtx_index][1],2)+pow(tang_pre[vtx_index][2],2));
00482         if (fabs(dotproduct) > length){
00483                 if (dotproduct > 0)
00484                         dotproduct = length;
00485                 else
00486                         dotproduct = -1.0*length;
00487         }
00488         theta1 = 180.0/3.1415926*acos(dotproduct/length);
00489         u = u3[vtx_index] + 0.5*(u4[vtx_index]-u3[vtx_index]);
00490         g_err = mk_core()->igeom_instance()->getEntUtoXYZ(edges[sorted_edge_list[next]].gEdgeHandle, u, pts1[0], pts1[1], pts1[2]);
00491         IBERRCHK(g_err, "Trouble get xyz coordinates from parametric coordinates on edge.");
00492         u = u3[vtx_index] + 0.501*(u4[vtx_index]-u3[vtx_index]);
00493         g_err = mk_core()->igeom_instance()->getEntUtoXYZ(edges[sorted_edge_list[next]].gEdgeHandle, u, pts2[0], pts2[1], pts2[2]);
00494         IBERRCHK(g_err, "Trouble get xyz coordinates from parametric coordinates on edge.");
00495         for (int j = 0; j < 3; j++)
00496                 vec[j] = pts2[j] - pts1[j];
00497         dotproduct = vec[0]*tang_next[vtx_index][0] + vec[1]*tang_next[vtx_index][1] + vec[2]*tang_next[vtx_index][2];
00498         length = sqrt(pow(vec[0],2)+pow(vec[1], 2)+pow(vec[2], 2))*sqrt(pow(tang_next[vtx_index][0],2)+pow(tang_next[vtx_index][1],2)+pow(tang_next[vtx_index][2],2));
00499         if (fabs(dotproduct) > length){
00500                 if (dotproduct > 0)
00501                         dotproduct = length;
00502                 else
00503                         dotproduct = -1.0*length;
00504         }
00505         theta2 = 180.0/3.1415926*acos(dotproduct/length);
00506         
00507         if ((fabs(theta1) > 0.5)&&(fabs(theta2) > 0.5))
00508                 return true;
00509         else
00510                 return false;
00511 
00512 }
00513 
00514 
00515 //reorganize the vertices and edges on the boudnaries of geometry
00516 void SubMapping::VerEdgOrganize(std::set<iBase_EntityHandle> edge_set, std::vector<iBase_EntityHandle> g_edge, iBase_EntityHandle surf)
00517 {
00518         int first_edge_index = 0;
00519         int first_index = edges[0].connect[1]->index;
00520         int second_index = edges[0].connect[0]->index;
00521         int start_index = edges[0].connect[0]->index;
00522 
00523         int test_sense = -10;
00524         iGeom::Error g_err = mk_core()->igeom_instance()->getEgVtxSense(edges[first_edge_index].gEdgeHandle, vertices[second_index].gVertexHandle, vertices[first_index].gVertexHandle, test_sense);
00525         IBERRCHK(g_err, "Trouble get the sense of edge with respect to two vertices.");
00526 
00527         int test_face_sense = -10;
00528         g_err = mk_core()->igeom_instance()->getEgFcSense(edges[first_edge_index].gEdgeHandle, surf, test_face_sense);
00529         IBERRCHK(g_err, "Trouble get the sense of edge with respect to two vertices.");
00530         std::cout << "edge sense = " << test_sense << "\tface sense = " << test_face_sense << "combined(multiply) = " << test_sense*test_face_sense << std::endl;       
00531 
00532         if ((test_sense*test_face_sense) < 0){
00533                 int tmp = second_index;
00534                 second_index = first_index;
00535                 first_index = tmp;
00536                 start_index = second_index;
00537         }
00538         
00539 
00540         sorted_vertex_list.push_back(start_index);
00541         
00542         sorted_edge_list.push_back(first_edge_index);
00543 
00544         while(start_index != first_index){
00545                 sorted_vertex_list.push_back(first_index);
00546 
00547                 vector<iBase_EntityHandle> adj_edges;
00548                 g_err = mk_core()->igeom_instance()->getEntAdj(vertices[first_index].gVertexHandle, iBase_EDGE, adj_edges);
00549                 IBERRCHK(g_err, "Trouble get the adjacent edges w.r.t a vertex.");
00550 
00551                 int index = -1;
00552                 for (unsigned int i = 0; i < adj_edges.size(); i++){
00553                         if ((adj_edges[i] != edges[first_edge_index].gEdgeHandle)&&(edge_set.find(adj_edges[i]) != edge_set.end())){
00554                                 g_err = mk_core()->igeom_instance()->getIntData(adj_edges[i], g_taghandle, index);
00555                                 IBERRCHK(g_err, "Trouble get the int data for geometric edge.");
00556                                 break;
00557                         }
00558                         else
00559                                 continue;
00560                 }
00561                 
00562                 //insert the new edge into sorted_edge_list
00563                 if (index > -1){                
00564                         
00565                         first_edge_index = index;
00566                         sorted_edge_list.push_back(first_edge_index);
00567                         second_index = first_index;
00568                         if (edges[first_edge_index].connect[0]->index == first_index)
00569                                 first_index = edges[first_edge_index].connect[1]->index;
00570                         else
00571                                 first_index = edges[first_edge_index].connect[0]->index;
00572                         
00573                 }               
00574                 else
00575                         break;  
00576         }
00577 }
00578 
00579 
00580 //classify the boundary edges as -I, +I, -J, +J
00581 void SubMapping::EdgeClassification()
00582 {
00583         edges_types.resize(sorted_edge_list.size());
00584         
00585         //Find a starting END vertex
00586         start_index = 0;
00587         for (; start_index < sorted_vertex_list.size(); start_index++)
00588                 if (vertices_types[sorted_vertex_list[start_index]]==END)
00589                         break;
00590 
00591         //setting up the initial direction for i-j coordinate system(positive i and positive j direction)
00592         int pre_edge = (start_index + sorted_edge_list.size() -1)%sorted_edge_list.size(), next_edge = start_index;
00593                 
00594         edges_types[sorted_edge_list[next_edge]] = POSI_J;
00595         edges_types[sorted_edge_list[pre_edge]] = NEG_I;
00596         EdgeTypes i_direction = NEG_I, j_direction = POSI_J;
00597         VertexTypes pre_vertex_type = END;
00598 
00599         int vertex_index = start_index;
00600         for (unsigned int i = 1; i < sorted_vertex_list.size()-1; i++){
00601                 pre_edge = (start_index + i + sorted_edge_list.size() -1)%sorted_edge_list.size(); //index in the sorted_edge_list
00602                 next_edge = (start_index + i)%sorted_edge_list.size();   //index in the sorted_edge_list
00603                 vertex_index = (i + start_index)%sorted_vertex_list.size();
00604 
00605                 switch(vertices_types[sorted_vertex_list[vertex_index]]){
00606                         case END://switch from i to j or from j to i
00607                                 if ((edges_types[sorted_edge_list[pre_edge]] == POSI_I)||(edges_types[sorted_edge_list[pre_edge]] == NEG_I)){
00608                                                 if (j_direction == POSI_J)
00609                                                         if (pre_vertex_type == END)
00610                                                                 edges_types[sorted_edge_list[next_edge]] = NEG_J;
00611                                                         else
00612                                                                 edges_types[sorted_edge_list[next_edge]] = POSI_J;
00613                                                 else
00614                                                         if (pre_vertex_type == END)
00615                                                                 edges_types[sorted_edge_list[next_edge]] = POSI_J;
00616                                                         else
00617                                                                 edges_types[sorted_edge_list[next_edge]] = NEG_J;
00618                                                 j_direction = edges_types[sorted_edge_list[next_edge]];                                 
00619                                 }
00620                                 else{//switch from j-direction to i-direction
00621                                         if (i_direction == POSI_I)
00622                                                 if (pre_vertex_type == END)
00623                                                         edges_types[sorted_edge_list[next_edge]] = NEG_I;
00624                                                 else
00625                                                         edges_types[sorted_edge_list[next_edge]] = POSI_I;
00626                                         else
00627                                                 if (pre_vertex_type == END)
00628                                                         edges_types[sorted_edge_list[next_edge]] = POSI_I;
00629                                                 else
00630                                                         edges_types[sorted_edge_list[next_edge]] = NEG_I;
00631                                         i_direction = edges_types[sorted_edge_list[next_edge]];
00632                                 }
00633                                 pre_vertex_type = END;                          
00634                                 break;
00635                         case SIDE://keep the same as previous
00636                                 edges_types[sorted_edge_list[next_edge]] = edges_types[sorted_edge_list[pre_edge]];
00637                                 break;
00638                 
00639                         case CORNER://switch from i to j or from j to i
00640                                 if ((edges_types[sorted_edge_list[pre_edge]] == POSI_I)||(edges_types[sorted_edge_list[pre_edge]] == NEG_I)){
00641                                         if (j_direction == POSI_J)
00642                                                 if (pre_vertex_type == END)
00643                                                         edges_types[sorted_edge_list[next_edge]] = POSI_J;
00644                                                 else
00645                                                         edges_types[sorted_edge_list[next_edge]] = NEG_J;
00646                                         else
00647                                                 if (pre_vertex_type == END)
00648                                                         edges_types[sorted_edge_list[next_edge]] = NEG_J;
00649                                                 else
00650                                                         edges_types[sorted_edge_list[next_edge]] = POSI_J;
00651                                         j_direction = edges_types[sorted_edge_list[next_edge]];                                         
00652                                 }
00653                                 else{
00654                                         if (i_direction == POSI_I)
00655                                                 if (pre_vertex_type == END)
00656                                                         edges_types[sorted_edge_list[next_edge]] = POSI_I;
00657                                                 else
00658                                                         edges_types[sorted_edge_list[next_edge]] = NEG_I;
00659                                         else
00660                                                 if (pre_vertex_type == END)
00661                                                         edges_types[sorted_edge_list[next_edge]] = NEG_I;
00662                                                 else
00663                                                         edges_types[sorted_edge_list[next_edge]] = POSI_I;
00664                                         i_direction = edges_types[sorted_edge_list[next_edge]]; 
00665                                 }
00666                                 pre_vertex_type = CORNER;                               
00667                                 break;
00668 
00669                         case REVERSAL://need to consider this point later on
00670                                 if (edges_types[sorted_edge_list[pre_edge]] == POSI_I)
00671                                         edges_types[sorted_edge_list[next_edge]] = NEG_I;
00672                                 else if (edges_types[sorted_edge_list[pre_edge]] == NEG_I)
00673                                         edges_types[sorted_edge_list[next_edge]] = POSI_I;
00674                                 else if (edges_types[sorted_edge_list[pre_edge]] == POSI_J)
00675                                         edges_types[sorted_edge_list[next_edge]] = NEG_J;
00676                                 else
00677                                         edges_types[sorted_edge_list[next_edge]] = POSI_J;
00678                                 break;
00679 
00680                         default://do nothing
00681                                 break;
00682                 }
00683         }       
00684 }
00685 
00686 //assign the i,j coordinates for each node on the boundary
00687 void SubMapping::EdgeDiscretization(ModelEnt *me)
00688 {
00689         edge_size.resize(sorted_edge_list.size());
00690         for (unsigned int i = 0; i < sorted_edge_list.size(); i++){
00691                 double measure;
00692                 iGeom::Error g_err = mk_core()->igeom_instance()->measure(&(edges[sorted_edge_list[i]].gEdgeHandle), 1, &measure);
00693                 IBERRCHK(g_err, "Trouble measure the boundary edges.");
00694                 edge_size[sorted_edge_list[i]] = int(measure/size_low_bound);
00695         }
00696 
00697         //linear programming to get # of line segments for each boundary edge
00698         LPSolveClass lp;
00699         //setup the model for linear programming
00700         vector<vector<double> > coeffs;
00701         vector<double> b(edges.size(), 1.0);
00702         
00703         //setup the objective function for linear programming
00704         lp.SetupObj(b, 0.0);
00705 
00706         //setup the equality constraint for linear programming
00707         coeffs.resize(2);
00708         coeffs[0].resize(edges.size());
00709         coeffs[1].resize(edges.size());
00710         b.clear();
00711         b.resize(2, 0.0);
00712         for (unsigned int i = 0; i < edges.size(); i++){
00713                 if ((edges_types[i] == POSI_I)||(edges_types[i] == NEG_I)){//positive i and negative i
00714                         if (edges_types[i] == POSI_I)
00715                                 coeffs[0][i] = 1.0;
00716                         else if (edges_types[i] == NEG_I)
00717                                 coeffs[0][i] = -1.0;
00718                         coeffs[1][i] = 0.0;
00719                 }
00720                 else{//positive j and negative j
00721                         if (edges_types[i] == POSI_J)
00722                                 coeffs[1][i] = 1.0;
00723                         else if (edges_types[i] == NEG_J)
00724                                 coeffs[1][i] = -1.0;
00725                         coeffs[0][i] = 0.0;
00726                 }
00727         }
00728         lp.SetupEqu(coeffs, b);
00729         
00730         //setup the constant constraint there is already mesh on the geometric edge
00731         vector<int> num_line_segments(edges.size());
00732         for (unsigned int i = 0; i < edges.size(); i++){
00733                 iBase_EntitySetHandle entityset;
00734                 iRel::Error r_err = mk_core()->irel_pair()->getEntSetRelation(edges[i].gEdgeHandle, 0, entityset);
00735                 IBERRCHK(r_err, "Trouble get the entity set for geometric edge.");
00736                 iMesh::Error m_err = mk_core()->imesh_instance()->getNumOfType(entityset, iBase_EDGE, num_line_segments[i]);
00737                 IBERRCHK(m_err, "Trouble get # of line segments on the edge mesh.");
00738 
00739                 if ((num_line_segments[i] < 1)&&(!m_err))
00740                         num_line_segments[i] = -1;
00741         }
00742         lp.SetupConst(num_line_segments);
00743 
00744         //setup the inequality constraint for linear programming
00745         coeffs.clear();
00746         coeffs.resize(edges.size(), vector<double>(edges.size()));
00747         for (unsigned int i = 0; i < edges.size(); i++)//diagonal matrix
00748                 coeffs[i][i] = -1.0;
00749         b.clear();
00750         b.resize(edges.size());
00751         int min_line_segment = 1.0e10;
00752         for (unsigned int i = 0; i < num_line_segments.size(); i++)
00753                 if ((num_line_segments[i] < min_line_segment)&&(num_line_segments[i]>0))
00754                         min_line_segment = num_line_segments[i];
00755 
00756         //preprocess the inequality constraints
00757         int sum_posi_I = 0, sum_neg_I = 0, sum_posi_J = 0, sum_neg_J = 0;
00758         vector<int> posi_i, neg_i, posi_j, neg_j;
00759         for (int i = 0; i < num_line_segments.size(); i++){
00760                 if (num_line_segments[i] > 0){
00761                         switch(edges_types[i]){
00762                                 case POSI_I:
00763                                         sum_posi_I += num_line_segments[i]; break;
00764                                 case NEG_I:
00765                                         sum_neg_I += num_line_segments[i]; break;
00766                                 case POSI_J:
00767                                         sum_posi_J += num_line_segments[i]; break;
00768                                 case NEG_J:
00769                                         sum_neg_J += num_line_segments[i]; break;
00770                                 default:
00771                                         break;
00772                         }
00773                 }
00774                 else{
00775                         switch(edges_types[i]){
00776                                 case POSI_I:
00777                                         posi_i.push_back(i); break;
00778                                 case NEG_I:
00779                                         neg_i.push_back(i); break;
00780                                 case POSI_J:
00781                                         posi_j.push_back(i); break;
00782                                 case NEG_J:
00783                                         neg_j.push_back(i); break;
00784                                 default:
00785                                         break;          
00786                         }
00787                 }
00788         }
00789         //preprocess the edge size
00790         int tmp_sum_POS_I = sum_posi_I, tmp_sum_POS_J = sum_posi_J, tmp_sum_NEG_I = sum_neg_I, tmp_sum_NEG_J = sum_neg_J;
00791         for (unsigned int i = 0; i < posi_i.size(); i++)
00792                 tmp_sum_POS_I += edge_size[posi_i[i]];
00793         for (unsigned int i = 0; i < posi_j.size(); i++)
00794                 tmp_sum_POS_J += edge_size[posi_j[i]];
00795         for (unsigned int i = 0; i < neg_i.size(); i++)
00796                 tmp_sum_NEG_I += edge_size[neg_i[i]];
00797         for (unsigned int i = 0; i < neg_j.size(); i++)
00798                 tmp_sum_NEG_J += edge_size[neg_j[i]];
00799         
00800         //process i-direction ----new
00801         if (tmp_sum_POS_I == tmp_sum_NEG_I)
00802         {}//do nothing, it is always true
00803         else if (tmp_sum_POS_I > tmp_sum_NEG_I){//adjust the NEG_I
00804                 if (neg_i.size() > 0){
00805                         for (unsigned int i = 0; i < neg_i.size(); i++)
00806                                 edge_size[neg_i[i]] = int(double(tmp_sum_POS_I-sum_neg_I)*double(edge_size[neg_i[i]])/double(tmp_sum_NEG_I-sum_neg_I));
00807                 }
00808                 else{//neg_i.size == 0
00809                         if (posi_i.size() == 0){
00810                                 std::cout << "Constraint check fails in i-direction_a\n";exit(1);
00811                         }
00812                         else{
00813                                 if (sum_posi_I >= tmp_sum_NEG_I){//here, tmp_sum_NEG_I = sum_neg_I
00814                                         std::cout << "Constraint check fails in i-direction_b\n";exit(1);
00815                                 }
00816                                 else{//sum_posi_I < tmp_sum_NEG_I, adjust the edge posi_i
00817                                         for (unsigned int i = 0; i < posi_i.size(); i++){
00818                                                 edge_size[posi_i[i]] = int(double(tmp_sum_NEG_I-sum_posi_I)*double(edge_size[posi_i[i]])/double(tmp_sum_POS_I-sum_posi_I));
00819                                         }
00820                                 }
00821                         }
00822                 }
00823         }
00824         else{//tmp_sum_POS_I < tmp_sum_NEG_I, //adjust the POSS_I
00825                 if (posi_i.size() > 0){
00826                         for (unsigned int i = 0; i < posi_i.size(); i++)
00827                                 edge_size[posi_i[i]] = int(double(tmp_sum_NEG_I-sum_posi_I)*double(edge_size[posi_i[i]])/double(tmp_sum_POS_I-sum_posi_I));
00828                 }
00829                 else{//posi_i.size() == 0
00830                         //further constraints, 
00831                         if (neg_i.size() == 0){
00832                                 std::cout << "Constraint check fails in i-direction_c\n";exit(1);
00833                         }
00834                         else{
00835                                 if (sum_neg_I >= tmp_sum_POS_I){//here tmp_sum_POS_I = sum_posi_I
00836                                         std::cout << "Constraint check fails in i-direction_d\n";exit(1);
00837                                 }
00838                                 else{//sum_neg_I < tmp_sum_POS_I, adjust the edge neg_i
00839                                         for (unsigned int i = 0; i < neg_i.size(); i++){
00840                                                 edge_size[neg_i[i]] = int(double(tmp_sum_POS_I-sum_neg_I)*double(edge_size[neg_i[i]])/double(tmp_sum_NEG_I-sum_neg_I));
00841                                         }
00842                                 }
00843                         }
00844                 }
00845         }
00846 
00847 
00848         
00849         //process j-direction ----new
00850         if (tmp_sum_POS_J == tmp_sum_NEG_J)
00851         {}//do nothing, it is always true
00852         else if (tmp_sum_POS_J > tmp_sum_NEG_J){//adjust the NEG_J
00853                 if (neg_j.size() > 0){
00854                         for (unsigned int i = 0; i < neg_j.size(); i++)
00855                                 edge_size[neg_j[i]] = int(double(tmp_sum_POS_J-sum_neg_J)*double(edge_size[neg_j[i]])/double(tmp_sum_NEG_J-sum_neg_J));
00856                 }
00857                 else{//neg_j.size == 0
00858                         if (posi_j.size() == 0){
00859                                 std::cout << "Constraint check fails in j-direction_a\n";exit(1);
00860                         }
00861                         else{
00862                                 if (sum_posi_J >= tmp_sum_NEG_J){//here, tmp_sum_NEG_J = sum_neg_J
00863                                         std::cout << "Constraint check fails in j-direction_b\n";exit(1);
00864                                 }
00865                                 else{//sum_posi_J < tmp_sum_NEG_J, adjust the edge posi_j
00866                                         for (unsigned int i = 0; i < posi_j.size(); i++){
00867                                                 edge_size[posi_j[i]] = int(double(tmp_sum_NEG_J-sum_posi_J)*double(edge_size[posi_j[i]])/double(tmp_sum_POS_J-sum_posi_J));
00868                                         }
00869                                 }
00870                         }
00871                 }
00872         }
00873         else{//tmp_sum_POS_J < tmp_sum_NEG_J, //adjust the POS_I
00874                 if (posi_j.size() > 0){
00875                         for (unsigned int i = 0; i < posi_j.size(); i++)
00876                                 edge_size[posi_j[i]] = int(double(tmp_sum_NEG_J-sum_posi_J)*double(edge_size[posi_j[i]])/double(tmp_sum_POS_J-sum_posi_J));
00877                 }
00878                 else{//posi_j.size() == 0
00879                         //further constraints, 
00880                         if (neg_j.size() == 0){
00881                                 std::cout << "Constraint check fails in j-direction_c\n";exit(1);
00882                         }
00883                         else{
00884                                 if (sum_neg_J >= tmp_sum_POS_J){//here tmp_sum_POS_J = sum_posi_J
00885                                         std::cout << "Constraint check fails in j-direction_d\n";exit(1);
00886                                 }
00887                                 else{//sum_neg_J < tmp_sum_POS_J, adjust the edge neg_j
00888                                         for (unsigned int i = 0; i < neg_j.size(); i++){
00889                                                 edge_size[neg_j[i]] = int(double(tmp_sum_POS_J-sum_neg_J)*double(edge_size[neg_j[i]])/double(tmp_sum_NEG_J-sum_neg_J));
00890                                         }
00891                                 }
00892                         }
00893                 }
00894         }
00895 
00896 
00897 
00898         
00899 
00900         //process i-direction
00901         if ((sum_posi_I == sum_neg_I)&&(sum_posi_I != 0)){
00902                 if ((posi_i.size() != 0)&&(neg_i.size() != 0))
00903                 {}//these constraints should work
00904                 else if ((posi_i.size() == 0)&&(neg_i.size() == 0))
00905                 {}//these constraints should work
00906                 else if ((posi_i.size() != 0)&&(neg_i.size() == 0))//this constraints won't work
00907                 {std::cout << "Constraint check fails in i-direction3\n";exit(1);}
00908                 else//these constraints won't work
00909                 {std::cout << "Constraint check fails in i-direction4\n";exit(1);}
00910         }
00911         else{
00912                 if (sum_posi_I > sum_neg_I){
00913                         if (sum_neg_I == 0){//low bound for NEG_I edges should be less than {sum_posi_I+edge_size[POSI_I]}
00914                                 //this work     
00915                         }
00916                         else{//sum_neg_I > 0
00917                                 if ((int(posi_i.size()) == 0)&&(int(neg_i.size()) == 0))//these constraints won't work
00918                                 {std::cout << "Constraint check fails in i-direction5\n";exit(1);}
00919                                 else if ((int(posi_i.size()) == 0)&&(int(neg_i.size()) > 0)){//sum_neg_I + edge_size[NEG_I] <= sum_posi_I
00920                                         //this works
00921                                 }
00922                                 else if ((int(posi_i.size()) > 0)&&(int(neg_i.size()) == 0))//these constraints won't work
00923                                 {std::cout << "Constraint check fails in i-direction6\n";exit(1);}
00924                                 else{ //((int(posi_i.size()) > 0)&&(int(neg_i.size()) > 0))
00925                                         //these constraints should work.
00926                                 }
00927                         }
00928                 }
00929                 else{//sum_posi_I < sum_neg_I
00930                         if (sum_posi_I == 0){//low bound for POSI_I edges should be less than {sum_neg_I+edge_size[NEG_I]}
00931                                 //this works    
00932                         }
00933                         else{//sum_posi_I > 0
00934                                 if ((int(posi_i.size()) == 0)&&(int(neg_i.size()) == 0))//these constraints won't work
00935                                 {std::cout << "Constraint check fails in i-direction7\n";exit(1);}
00936                                 else if ((int(neg_i.size()) == 0)&&(int(posi_i.size()) > 0)){//sum_posi_I + edge_size[POSI_I] <= sum_neg_I
00937                                         //this works
00938                                 }
00939                                 else if ((int(neg_i.size()) > 0)&&(int(posi_i.size()) == 0))//these constraints won't work
00940                                 {std::cout << "Constraint check fails in i-direction8\n";exit(1);}
00941                                 else{ //((int(posi_i.size()) > 0)&&(int(neg_i.size()) > 0))
00942                                         //these constraints should work.
00943                                 }
00944                         }
00945                 }
00946         }
00947         //process j direction
00948         if ((sum_posi_J == sum_neg_J)&&(sum_posi_J != 0)){
00949                 if ((posi_j.size() != 0)&&(neg_j.size() != 0))//these constraints should work
00950                 {}
00951                 else if ((posi_j.size() == 0)&&(neg_j.size() == 0))
00952                 {}//these constraints should work
00953                 else if ((posi_j.size() != 0)&&(neg_j.size() == 0))//this constraints won't work
00954                 {std::cout << "Constraint check fails in j-direction3\n";exit(1);}
00955                 else//these constraints won't work
00956                 {std::cout << "Constraint check fails in j-direction4\n";exit(1);}
00957         }
00958         else{//sum_posi_J != sum_neg_J
00959                 if (sum_posi_J > sum_neg_J){
00960                         if (sum_neg_J == 0){//low bound for NEG_J edges should be less than {sum_posi_J+edge_size[POSI_J]}
00961                                 //this works    
00962                         }
00963                         else{//sum_neg_J > 0
00964                                 if ((int(posi_j.size()) == 0)&&(int(neg_j.size()) == 0))//these constraints won't work
00965                                 {std::cout << "Constraint check fails in j-direction5\n";exit(1);}
00966                                 else if ((int(posi_j.size()) == 0)&&(int(neg_j.size()) > 0)){//sum_neg_J + edge_size[NEG_J] <= sum_posi_J
00967                                         //this works
00968                                 }
00969                                 else if ((int(posi_j.size()) > 0)&&(int(neg_j.size()) == 0))//these constraints won't work
00970                                 {std::cout << "Constraint check fails in j-direction6\n";exit(1);}
00971                                 else{ //((int(posi_j.size()) > 0)&&(int(neg_j.size()) > 0))
00972                                         //these constraints should work.
00973                                 }
00974                         }
00975                 }
00976                 else{//sum_posi_J < sum_neg_J
00977                         if (sum_posi_J == 0){//low bound for POSI_J edges should be less than {sum_neg_J+edge_size[NEG_J]}
00978                                 //this works            
00979                         }
00980                         else{//sum_posi_J > 0
00981                                 if ((int(posi_j.size()) == 0)&&(int(neg_j.size()) == 0))//these constraints won't work
00982                                 {std::cout << "Constraint check fails in j-direction7\n";exit(1);}
00983                                 else if ((int(neg_j.size()) == 0)&&(int(posi_j.size()) > 0)){//sum_posi_J + edge_size[POSI_J] <= sum_neg_J
00984                                         //this works
00985                                 }
00986                                 else if ((int(neg_j.size()) > 0)&&(int(posi_j.size()) == 0))//these constraints won't work
00987                                 {std::cout << "Constraint check fails in j-direction8\n";exit(1);}
00988                                 else{ //((int(posi_j.size()) > 0)&&(int(neg_j.size()) > 0))
00989                                         //these constraints should work.
00990                                 }
00991                         }
00992                 }
00993         }
00994         
00995         for (unsigned int i = 0; i < edges.size(); i++)
00996                 if (min_line_segment < 1.0e9){
00997                         if ((num_line_segments[i] < edge_size[i])&&(num_line_segments[i] > 0)){//avoid the conflicts between const constraints and inequality constraints;
00998                                 b[i] = -1.0*num_line_segments[i];
00999                         }
01000                         else{
01001                                 //if (min_line_segment < edge_size[i])
01002                                 //      b[i] = -1.0*min_line_segment;
01003                                 //else
01004                                 b[i] = -1.0*edge_size[i];
01005                         }
01006                 }
01007                 else
01008                         b[i] = -1.0*edge_size[i];
01009         
01010         lp.SetupInEqu(coeffs, b);
01011 
01012         lp.Execute();
01013         lp.GetVariables(edge_size);
01014 
01015         //do the vertex mesher
01016         //get the mesh node on the geometric vertex
01017         for (unsigned int i = 0; i < sorted_vertex_list.size(); i++){
01018                 iBase_EntitySetHandle entityset;
01019                 iRel::Error r_err = mk_core()->irel_pair()->getEntSetRelation(vertices[sorted_vertex_list[i]].gVertexHandle, 0, entityset);
01020                 IBERRCHK(r_err, "Trouble get the entity set for edge from geometric vertex.");
01021 
01022                 vector<iBase_EntityHandle> points;
01023                 points.clear();
01024                 iMesh::Error m_err = mk_core()->imesh_instance()->getEntities(entityset, iBase_VERTEX, iMesh_POINT, points);
01025                 IBERRCHK(m_err, "Trouble get the node entities from mesh entity sets.");
01026 
01027                 if (points.size()==0){
01028                         points.resize(1);
01029                         //there is no mesh nodes on the geometric vertices, it need to create a mesh node on the geometric vertex
01030                         m_err = mk_core()->imesh_instance()->createVtx(vertices[sorted_vertex_list[i]].xyz[0], vertices[sorted_vertex_list[i]].xyz[1], vertices[sorted_vertex_list[i]].xyz[2], points[0]);
01031                         IBERRCHK(m_err, "Trouble create the mesh nodes on the geometric vertex.");
01032                         m_err = mk_core()->imesh_instance()->addEntToSet(points[0], entityset);
01033                         IBERRCHK(m_err, "Trouble add the mesh node entity to entity set.");
01034                 }
01035         }
01036 
01037         //do the edge mesher
01038         MEntVector curves;
01039         me->get_adjacencies(1, curves);
01040         assert(curves.size()==edges.size());
01041 
01042         //call the edge mesher to discretize the boundary edges
01043         for (unsigned int i = 0; i < edges.size(); i++){
01044                 //initial size functon for edges, get the number of edges and assign it to the edge
01045 
01046                 //do the edge mesher 
01047                 SizingFunction esize(mk_core(), edge_size[i], -1);
01048                 me->sizing_function_index(esize.core_index());
01049 
01050                 //detect the edge on the surface
01051                 MEntVector edge_curve;
01052                 edge_curve.resize(1);
01053                 
01054                 for (unsigned int j = 0; j < curves.size(); j++){
01055                         int index_id;
01056                         iGeom::Error g_err = mk_core()->igeom_instance()->getIntData(curves[j]->geom_handle(), g_taghandle, index_id);
01057                         IBERRCHK(g_err, "Trouble get the int data for geometric edges.");
01058 
01059                         if (index_id == (int)i){
01060                                 edge_curve[0] = curves[j];
01061                                 break;
01062                         }
01063                 }
01064                 
01065                 //do the edge mesher
01066                 EdgeMesher *em = (EdgeMesher*) mk_core()->construct_meshop("EdgeMesher", edge_curve);
01067                 //mk_core()->setup_and_execute();
01068                 em->setup_this();
01069                 em->execute_this();     
01071         }
01072 
01073         //extract the mesh info
01074         iMesh::Error m_err = mk_core()->imesh_instance()->getTagHandle("MeshSubMapping", m_taghandle);
01075         if (m_err){
01076                 m_err = mk_core()->imesh_instance()->createTag("MeshSubMapping", 1, iBase_INTEGER, m_taghandle);        
01077                 IBERRCHK(m_err, "Trouble create a taghandle.");
01078         }
01079 
01080         int index = 0;  
01081         for (unsigned int i = 0; i < sorted_vertex_list.size(); i++){
01082                 int next_edge = sorted_edge_list[i%sorted_edge_list.size()];
01083                 int pre_edge = sorted_edge_list[(i+sorted_edge_list.size()-1)%sorted_edge_list.size()];
01084                 int next_vertex = sorted_vertex_list[(i+1)%sorted_vertex_list.size()];
01085                 iBase_EntitySetHandle entityset;
01086 
01087                 //get the mesh node on the geometric vertex
01088                 iRel::Error r_err = mk_core()->irel_pair()->getEntSetRelation(vertices[sorted_vertex_list[i]].gVertexHandle, 0, entityset);
01089                 IBERRCHK(r_err, "Trouble get the entity set for edge from geometric vertex.");
01090 
01091                 vector<iBase_EntityHandle> points;
01092                 points.clear();
01093                 m_err = mk_core()->imesh_instance()->getEntities(entityset, iBase_VERTEX, iMesh_POINT, points);
01094                 IBERRCHK(m_err, "Trouble get the node entities from mesh entity sets.");
01095 
01096                 assert(points.size()==1);
01097                 index++;
01098                 nodes.resize(index);
01099                 nodes[index-1].index = index-1;
01100                 nodes[index-1].onBoundary = false;
01101                 nodes[index-1].onCorner = true;
01102                 geom_mesh_node[sorted_vertex_list[i]] = index -1;
01103                 mesh_geom_vertex[index-1] = sorted_vertex_list[i];
01104                 nodes[index-1].gVertexHandle = points[0];
01105                 m_err = mk_core()->imesh_instance()->setIntData(nodes[index-1].gVertexHandle, m_taghandle, index-1);
01106                 IBERRCHK(m_err, "Trouble set the int data for mesh nodes.");
01107                 for (int j = 0; j < 3; j++)
01108                         nodes[index-1].xyz[j] = vertices[sorted_vertex_list[i]].xyz[j];
01109 
01110                 //assign the i-j coordinates for boundary nodes
01111                 if (i == 0){
01112                         nodes[index-1].uv[0] = 0.0;
01113                         nodes[index-1].uv[1] = 0.0;
01114                 }
01115                 else{
01116                         AssignIJCoords(nodes[index-1].uv[0], nodes[index-1].uv[1], edges_types[pre_edge], index -1);
01117                 }
01118 
01119                 //get the mesh nodes on the geometric edges
01120                 r_err = mk_core()->irel_pair()->getEntSetRelation(edges[next_edge].gEdgeHandle, 0, entityset);
01121                 IBERRCHK(r_err, "Trouble get the entity set for edge from geometric vertex.");
01122 
01123                 points.clear();
01124                 m_err = mk_core()->imesh_instance()->getEntities(entityset, iBase_VERTEX, iMesh_POINT, points);
01125                 IBERRCHK(m_err, "Trouble get the node entities from mesh entity sets.");
01126 
01127                 int sense = 0;
01128                 iGeom::Error g_err = mk_core()->igeom_instance()->getEgVtxSense(edges[next_edge].gEdgeHandle, vertices[sorted_vertex_list[i]].gVertexHandle, vertices[next_vertex].gVertexHandle, sense);
01129                 IBERRCHK(g_err, "Trouble get the edge sense with respect to two vertices.");
01130 
01131                 unsigned int j;
01132                 if (sense < 0) 
01133                         j = points.size()-1;
01134                 else
01135                         j = 0; 
01136                 for (; ((j < points.size())&&(j>=0));){                                         
01137                         index++;
01138                         nodes.resize(index);
01139                         nodes[index-1].index = index-1;
01140                         nodes[index-1].gVertexHandle = points[j];
01141                         nodes[index-1].onBoundary = true;
01142                         nodes[index-1].onCorner = false;
01143                         m_err = mk_core()->imesh_instance()->setIntData(nodes[index-1].gVertexHandle, m_taghandle, index-1);
01144                         IBERRCHK(m_err, "Trouble set the int data for mesh nodes.");
01145                         m_err = mk_core()->imesh_instance()->getVtxCoord(nodes[index-1].gVertexHandle, nodes[index-1].xyz[0], nodes[index-1].xyz[1], nodes[index-1].xyz[2]);
01146                         IBERRCHK(m_err, "Trouble get the coordinates from mesh nodes.");
01147                         AssignIJCoords(nodes[index-1].uv[0], nodes[index-1].uv[1], edges_types[next_edge], index -1);
01148                         if (sense < 0) j--;
01149                         else j++;
01150                 }
01151         }
01152 }
01153 
01154 //Assign the i j coordinates for boundary nodes
01155 void SubMapping::AssignIJCoords(double &u, double &v, EdgeTypes type, int index){
01156         switch(type){
01157                 case POSI_I:
01158                         nodes[index].uv[0] = nodes[index-1].uv[0]+1;
01159                         nodes[index].uv[1] = nodes[index-1].uv[1];
01160                         break;
01161                 case NEG_I:
01162                         nodes[index].uv[0] = nodes[index-1].uv[0]-1;
01163                         nodes[index].uv[1] = nodes[index-1].uv[1];
01164                         break;
01165                 case POSI_J:
01166                         nodes[index].uv[0] = nodes[index-1].uv[0];
01167                         nodes[index].uv[1] = nodes[index-1].uv[1]+1;
01168                         break;
01169                 case NEG_J:
01170                         nodes[index].uv[0] = nodes[index-1].uv[0];
01171                         nodes[index].uv[1] = nodes[index-1].uv[1]-1;
01172                         break;
01173                 default:
01174                         break;
01175         }       
01176 }
01177 
01178 //subdivide the surface 
01179 void SubMapping::InteriorNodeInterpolation(ModelEnt *me)
01180 {
01181         iBase_EntitySetHandle entityset;
01182         iRel::Error r_err = mk_core()->irel_pair()->getEntSetRelation(me->geom_handle(), 0, entityset);
01183         IBERRCHK(r_err, "Trouble get the entity set for edge from geometric vertex.");
01184 
01185         int i_max = -1.0e5, i_min = 1.0e5, j_max = -1.0e5, j_min = 1.0e5;
01186         //calculate the bounding i-j: i_max, i_min, j_max, j_min, these, these extreme points should be on the geometric vertices
01187         for (unsigned int i = 0; i < sorted_vertex_list.size(); i++){
01188                 int index = geom_mesh_node[sorted_vertex_list[i]];
01189                 if (nodes[index].uv[0] < i_min)
01190                         i_min = int(nodes[index].uv[0]);
01191                 if (nodes[index].uv[0] > i_max)
01192                         i_max = int(nodes[index].uv[0]);
01193                 if (nodes[index].uv[1] < j_min)
01194                         j_min = int(nodes[index].uv[1]);
01195                 if (nodes[index].uv[1] > j_max)
01196                         j_max = int(nodes[index].uv[1]);
01197         }
01198         vector<Vertex> boundary_ij;
01199         for (unsigned int i = 0; i < nodes.size(); i++)
01200                 if (nodes[i].onCorner){
01201                         boundary_ij.resize(boundary_ij.size()+1);
01202                         boundary_ij[boundary_ij.size()-1] = nodes[i];
01203         }
01204 
01205         //i_min, i_max, j_min and j_max are the boundaries
01206         int index = nodes.size();
01207         int num_node_boundary = index;
01208         for (int i = i_min+1; i < i_max; i++){
01209                 for (int j = j_min+1; j < j_max; j++){
01210                         double j_coord[2] = {j_min, j_max}, i_coord[2] = {i_min, i_max};
01211                         int j_index[2] = {-1, -1}, i_index[2] = {-1, -1};
01212                         for (int k = 0; k < num_node_boundary; k++){
01213                                 if ((int(nodes[k].uv[0]) == i)&&(int(nodes[k].uv[1]) == j)){
01214                                         //this point is a boundary node, i_index and j_index are always -1
01215                                         for (int m = 0; m < 2; m++){
01216                                                 i_index[m] = -1;
01217                                                 j_index[m] = -1;
01218                                         }
01219                                         break;
01220                                 }
01221                                 else if ((int(nodes[k].uv[0]) == i)&&(int(nodes[k].uv[1]) != j)){
01222                                         //this point has vertical nodes
01223                                         if ((nodes[k].uv[1] > j)&&(nodes[k].uv[1] <= j_coord[1])){//top point
01224                                                 j_index[1] = k;
01225                                                 j_coord[1] = nodes[k].uv[1];
01226                                         }
01227                                         else if ((nodes[k].uv[1] < j)&&(nodes[k].uv[1] >= j_coord[0])){//bottom point
01228                                                 j_index[0] = k;
01229                                                 j_coord[0] = nodes[k].uv[1];
01230                                         }
01231                                         else
01232                                                 continue;
01233                                 }
01234                                 else if ((int(nodes[k].uv[0]) != i)&&(int(nodes[k].uv[1]) == j)){
01235                                         //this point has horizontal nodes
01236                                         if ((nodes[k].uv[0] > i)&&(nodes[k].uv[0] <= i_coord[1])){
01237                                                 i_index[1] = k;
01238                                                 i_coord[1] = nodes[k].uv[0];
01239                                         }
01240                                         else if ((nodes[k].uv[0] < i)&&(nodes[k].uv[0] >= i_coord[0])){
01241                                                 i_index[0] = k;
01242                                                 i_coord[0] = nodes[k].uv[0];
01243                                         }
01244                                         else
01245                                                 continue;
01246                                 }
01247                                 else//this point is not interesting for me
01248                                         continue;
01249                         }
01250                         if ((i_index[0]!= -1)&&(i_index[1] != -1)&&(j_index[0] != -1)&&(j_index[1] != -1)){
01251                                 //interpolate the interior point
01252                                 double xyz[3];
01253                                 double weight[2] = {fabs(i-nodes[i_index[1]].uv[0]), fabs(j-nodes[j_index[1]].uv[1])};
01254                                 if (fabs(j-nodes[j_index[0]].uv[1]) < fabs(j-nodes[j_index[1]].uv[1]))
01255                                         weight[1] = double(fabs(j-nodes[j_index[0]].uv[1]));
01256                                 if (fabs(i-nodes[i_index[0]].uv[0]) < fabs(i-nodes[i_index[1]].uv[0]))
01257                                         weight[0] = fabs(i-nodes[i_index[0]].uv[0]);
01258                                 double total = weight[0]+weight[1];                             
01259                                 double i_weight = weight[1]/total, j_weight = weight[0]/total;
01260                                 
01261                                 for (int k = 0; k < 3; k++){
01262                                         double firstpart = nodes[j_index[1]].xyz[k]*fabs(j-nodes[j_index[0]].uv[1])/fabs(nodes[j_index[1]].uv[1]-nodes[j_index[0]].uv[1]);
01263                                         double secondpart = nodes[j_index[0]].xyz[k]*fabs(nodes[j_index[1]].uv[1]-j)/fabs(nodes[j_index[1]].uv[1]-nodes[j_index[0]].uv[1]);
01264                                         double j_value = firstpart + secondpart;
01265                                         double thirdpart = nodes[i_index[1]].xyz[k]*fabs(i-nodes[i_index[0]].uv[0])/fabs(nodes[i_index[1]].uv[0]-nodes[i_index[0]].uv[0]);
01266                                         double fourthpart = nodes[i_index[0]].xyz[k]*fabs(nodes[i_index[1]].uv[0]-i)/fabs(nodes[i_index[1]].uv[0]-nodes[i_index[0]].uv[0]);
01267                                         double i_value = thirdpart + fourthpart;
01268                                         xyz[k] = j_weight*firstpart+j_weight*secondpart+i_weight*thirdpart+i_weight*fourthpart;
01269                                 }
01270                                 iGeom::Error g_err = mk_core()->igeom_instance()->getEntClosestPt(me->geom_handle(), xyz[0], xyz[1], xyz[2], xyz[0], xyz[1], xyz[2]);
01271                                 IBERRCHK(g_err, "Trouble get the closest point on the surface.");
01272                                 if (isOutSideSurf(boundary_ij, i, j)){
01273                                 //if ((num_i[0]%2==1)||(num_i[1]%2==1)||(num_j[0]%2==1)||(num_j[1]%2==1))                       
01274                                         //create the vertex entity
01275                                         index++;
01276                                         nodes.resize(index);
01277                                         nodes[index-1].index = index - 1;
01278                                         for (int k = 0; k < 3; k++)
01279                                                 nodes[index-1].xyz[k] = xyz[k];
01280                                         nodes[index-1].onBoundary = false;
01281                                         nodes[index-1].onCorner = false;
01282                                         nodes[index-1].uv[0] = i;
01283                                         nodes[index-1].uv[1] = j;       
01284                                 }
01285                         }
01286                 }
01287         }
01288         
01289 
01290         std::cout << "number of boundary nodes is " << num_node_boundary << std::endl;
01291         //create the interior node entities
01292         vector<iBase_EntityHandle> m_nodes;
01293         m_nodes.resize(nodes.size()-num_node_boundary);
01294         for (unsigned int i = num_node_boundary; i < nodes.size(); i++){
01295                 iMesh::Error m_err = mk_core()->imesh_instance()->createVtx(nodes[i].xyz[0], nodes[i].xyz[1], nodes[i].xyz[2], m_nodes[i-num_node_boundary]);
01296                 IBERRCHK(m_err, "Trouble create the node entity.");
01297                 nodes[i].gVertexHandle = m_nodes[i-num_node_boundary];
01298                 m_err = mk_core()->imesh_instance()->addEntToSet(m_nodes[i-num_node_boundary], entityset);
01299                 IBERRCHK(m_err, "Trouble add the mesh node into the entityset.");
01300                 m_err = mk_core()->imesh_instance()->setIntData(m_nodes[i-num_node_boundary], m_taghandle, i);
01301                 IBERRCHK(m_err, "Trouble set the int data.");
01302         }
01303         //ok, we are done with the node entities
01304         
01305         //trick here, create an array of int data to quickly locate the node index
01306         vector<int> node_index((j_max-j_min+1)*(i_max-i_min+1),-1);//here, invalid position is -1
01307         for (unsigned int i = 0; i < nodes.size(); i++){
01308                 node_index[(nodes[i].uv[1]-j_min)*(i_max-i_min+1)+nodes[i].uv[0]-i_min] = i;
01309         }
01310 
01311         int face_sense;
01312         iGeom::Error g_err = mk_core()->igeom_instance()->getEgFcSense(edges[sorted_edge_list[0]].gEdgeHandle, me->geom_handle(), face_sense);
01313         IBERRCHK(g_err, "Trouble get the edge sense with respect to a geometrical face.");
01314 
01315         int edge_sense;
01316         g_err = mk_core()->igeom_instance()->getEgVtxSense(edges[sorted_edge_list[0]].gEdgeHandle, vertices[sorted_vertex_list[0]].gVertexHandle, vertices[sorted_vertex_list[1]].gVertexHandle, edge_sense);
01317         IBERRCHK(g_err, "Trouble get the sense of edge with respect to two vertices.");
01318 
01319         int test_sense = face_sense * edge_sense;
01320 
01321         index = 0;
01322         //create the quadrilateral elements on the surface
01323         for (int j = j_min; j < j_max; j++){
01324                 for (int i = i_min; i < i_max; i++){
01325                         
01326                         if (node_index[(j-j_min)*(i_max-i_min+1)+i-i_min] != -1){//there exists such a mesh node
01327                                 //check the other three vertices whether they exist or not
01328                                 if ((j+1)>j_max)
01329                                         continue;
01330                                 if ((i+1)>i_max)
01331                                         continue;
01332                                 if ((node_index[(j+1-j_min)*(i_max-i_min+1)+i-i_min] != -1)&&(node_index[(j-j_min)*(i_max-i_min+1)+i+1-i_min] != -1)&&(node_index[(j+1-j_min)*(i_max-i_min+1)+i+1-i_min] != -1)){
01333                                         index++;
01334                                         quads.resize(index);
01335                                         quads[index-1].connect.resize(4);
01336                                         if (test_sense < 0){
01337                                                 quads[index-1].connect[0] = &nodes[node_index[(j-j_min)*(i_max-i_min+1)+i-i_min]];
01338                                                 quads[index-1].connect[3] = &nodes[node_index[(j+1-j_min)*(i_max-i_min+1)+i-i_min]];
01339                                                 quads[index-1].connect[2] = &nodes[node_index[(j+1-j_min)*(i_max-i_min+1)+i+1-i_min]];
01340                                                 quads[index-1].connect[1] = &nodes[node_index[(j-j_min)*(i_max-i_min+1)+i+1-i_min]];
01341                                         }
01342                                         else{
01343                                                 quads[index-1].connect[0] = &nodes[node_index[(j-j_min)*(i_max-i_min+1)+i-i_min]];
01344                                                 quads[index-1].connect[1] = &nodes[node_index[(j+1-j_min)*(i_max-i_min+1)+i-i_min]];
01345                                                 quads[index-1].connect[2] = &nodes[node_index[(j+1-j_min)*(i_max-i_min+1)+i+1-i_min]];
01346                                                 quads[index-1].connect[3] = &nodes[node_index[(j-j_min)*(i_max-i_min+1)+i+1-i_min]];
01347                                         }                                       
01348                                 }
01349                         }
01350                 }
01351         }
01352         std::cout << "test quad size = " << quads.size() << std::endl;
01353 
01354         vector<iBase_EntityHandle> m_faces;
01355         m_faces.resize(quads.size());
01356         for (unsigned int i = 0; i < quads.size(); i++){
01357                 vector<iBase_EntityHandle> n_nodes;
01358                 n_nodes.resize(4);
01359                 for (int k = 0; k < 4; k++)
01360                         n_nodes[k] = nodes[quads[i].connect[k]->index].gVertexHandle;
01361 
01362                 iMesh::Error m_err = mk_core()->imesh_instance()->createEnt(iMesh_QUADRILATERAL, &n_nodes[0], 4, m_faces[i]);
01363                 IBERRCHK(m_err, "Trouble create the quadrilateral element.");
01364                 quads[i].gFaceHandle = m_faces[i];
01365                 m_err = mk_core()->imesh_instance()->addEntToSet(m_faces[i], entityset);
01366                 IBERRCHK(m_err, "Trouble add the mesh node into the entityset.");
01367                 m_err = mk_core()->imesh_instance()->setIntData(m_faces[i], m_taghandle, i);
01368                 IBERRCHK(m_err, "Trouble set the int data.");
01369 
01370                 n_nodes.clear();
01371         }
01372 
01373         std::cout << "Quad size = " << quads.size() << std::endl;       
01374         std::cout << "imax = " << i_max << "\timin = " << i_min << "\tjmax = " << j_max << "\tjmin = " << j_min << std::endl;
01375 }
01376 
01377 //smooth the structured mesh
01378 void SubMapping::MeshSmoothing(ModelEnt *ent)
01379 {
01380         std::vector<std::set<int> > AdjElements;
01381         std::vector<std::vector<int> > Quads(quads.size(), vector<int>(4));
01382         std::vector<std::vector<double> > coords(nodes.size(), vector<double>(3));
01383         std::vector<bool> isBnd(nodes.size(), false);
01384         std::vector<std::vector<int> > connect(nodes.size(), std::vector<int>(9));
01385         
01386         AdjElements.resize(nodes.size());
01387         for (unsigned int i = 0; i < quads.size(); i++){
01388             std::vector<iBase_EntityHandle> adjEnts;
01389             adjEnts.clear();
01390             iMesh::Error m_err = mk_core()->imesh_instance()->getEntAdj(quads[i].gFaceHandle, iBase_VERTEX, adjEnts);
01391             IBERRCHK(m_err, "Trouble get the adjacent nodes wrt a quad.");
01392 
01393             assert(adjEnts.size()==4);
01394 
01395             for (unsigned int j = 0; j < adjEnts.size(); j++){
01396                         int index_id = -1;
01397                         m_err = mk_core()->imesh_instance()->getIntData(adjEnts[j], m_taghandle, index_id);
01398                         IBERRCHK(m_err, "Trouble get int data for nodes.");
01399                         Quads[i][j] = index_id;
01400             }
01401         }
01402 
01403         for (unsigned int i = 0; i < nodes.size(); i++){
01404                 if (nodes[i].onBoundary || nodes[i].onCorner)
01405                         isBnd[i] = true;
01406                 for (int j = 0; j < 3; j++)
01407                         coords[i][j] = nodes[i].xyz[j];
01408                 if (!isBnd[i]){
01409                         vector<iBase_EntityHandle> adjEnts;
01410                         iMesh::Error m_err = mk_core()->imesh_instance()->getEntAdj(nodes[i].gVertexHandle, iBase_FACE, adjEnts);
01411                         IBERRCHK(m_err, "Trouble get the adjacent quads wrt a node.");
01412                         for (unsigned int j = 0; j < adjEnts.size(); j++){
01413                                 int index_id = -1;
01414                                 m_err = mk_core()->imesh_instance()->getIntData(adjEnts[j], m_taghandle, index_id);
01415                                 IBERRCHK(m_err, "Trouble get int data for quads.");
01416                                 AdjElements[i].insert(index_id);
01417                         }
01418 
01419                         //process the connect info
01420                         //there are 4 adjacent quadrilateral elements around node i     
01421                         assert(AdjElements[i].size() == 4);
01422                         std::set<int>::iterator it = AdjElements[i].begin();
01423                         int st_index[4];
01424                         //process 4 quad elements
01425                         int j = -1;
01426                         for (; it != AdjElements[i].end(); it++){
01427                                 j++;                            
01428                                 if (int(i) == Quads[*it][0])
01429                                         st_index[j] = 0;
01430                                 else if (int(i) == Quads[*it][1])
01431                                         st_index[j] = 1;
01432                                 else if (int(i) == Quads[*it][2])
01433                                         st_index[j] = 2;
01434                                 else
01435                                         st_index[j] = 3;
01436                         }
01437                         it = AdjElements[i].begin();
01438                         connect[i][2] = Quads[*it][(st_index[0]+3)%4];
01439                         connect[i][8] = Quads[*it][(st_index[0]+1)%4];
01440                         connect[i][1] = Quads[*it][(st_index[0]+2)%4];
01441                         //finish processing the quad 1
01442                         std::set<int>::iterator it1 = AdjElements[i].begin();
01443                         it1++;
01444                         for (j = 1; j < 4; j++, it1++){
01445                                 if (connect[i][8] == Quads[*it1][(st_index[j]+1)%4]){
01446                                         connect[i][7] = Quads[*it1][(st_index[j]+2)%4];
01447                                         connect[i][6] = Quads[*it1][(st_index[j]+3)%4];
01448                                         break;
01449                                 }
01450                                 else if (connect[i][8] == Quads[*it1][(st_index[j]+3)%4]){
01451                                         connect[i][7] = Quads[*it1][(st_index[j]+2)%4];
01452                                         connect[i][6] = Quads[*it1][(st_index[j]+1)%4];                                 
01453                                         break;
01454                                 }
01455                                 else
01456                                         continue;
01457                         }
01458                         //finish processing the quad 2
01459                         std::set<int>::iterator it2 = AdjElements[i].begin();
01460                         it2++;
01461                         for (j=1; it2 != AdjElements[i].end(); it2++,j++){
01462                                 if (connect[i][2] == Quads[*it2][(st_index[j]+1)%4]){
01463                                         connect[i][3] = Quads[*it2][(st_index[j]+2)%4];
01464                                         connect[i][4] = Quads[*it2][(st_index[j]+3)%4];
01465                                         break;
01466                                 }
01467                                 else if (connect[i][2] == Quads[*it2][(st_index[j]+3)%4]){
01468                                         connect[i][3] = Quads[*it2][(st_index[j]+2)%4];
01469                                         connect[i][4] = Quads[*it2][(st_index[j]+1)%4];                                 
01470                                         break;
01471                                 }
01472                                 else
01473                                         continue;
01474                         }
01475                         //finish processing the quad 4;
01476                         std::set<int>::iterator it3 = AdjElements[i].begin();
01477                         it3++;
01478                         for (j=1; it3 != AdjElements[i].end(); it3++,j++){
01479                                 if ((it3 != it1)&&(it3 != it2)){
01480                                         connect[i][5] = Quads[*it2][(st_index[j]+2)%4]; 
01481                                 }
01482                                 else
01483                                         continue;
01484                         }       
01485                 }
01486         }
01487 
01488         mk_core()->save_mesh("BeforeEquipotential.vtk");
01489         
01490         EquipotentialSmooth smoother;   
01491                 
01492         smoother.SetupData(AdjElements, Quads, coords, isBnd, connect);
01493         smoother.Execute();
01494         
01495         std::vector<std::vector<double> > coors(nodes.size(), vector<double>(3));
01496         smoother.GetCoords(coors);
01497 
01498         //update the new position for nodes
01499         for (unsigned int i = 0; i < nodes.size(); i++){
01500                 if (!isBnd[i]){
01501                         double tmp_coord[3] = {coords[i][0], coords[i][1], coords[i][2]};
01502                         
01503                         iGeom::Error g_err = mk_core()->igeom_instance()->getEntClosestPt(ent->geom_handle(), coords[i][0], coords[i][1], coords[i][2], tmp_coord[0], tmp_coord[1], tmp_coord[2]);
01504                         IBERRCHK(g_err, "Trouble get a closest point on the linking surface.");
01505                         iMesh::Error m_err = mk_core()->imesh_instance()->setVtxCoord(nodes[i].gVertexHandle, tmp_coord[0], tmp_coord[1], tmp_coord[2]);
01506                         IBERRCHK(m_err, "Trouble set the new coordinates for nodes.");
01507                 }
01508         }
01509         
01510 
01511 #ifdef HAVE_MESQUITE
01512         iBase_EntitySetHandle entityset;
01513         iRel::Error r_err = mk_core()->irel_pair()->getEntSetRelation(ent->geom_handle(), 0, entityset);
01514         IBERRCHK(r_err, "Trouble get the entity set for edge from geometric vertex.");
01515         MeshImprove shapesmooth(mk_core(), true, true, false, false, false);
01516     shapesmooth.SurfMeshImprove(ent->geom_handle(), entityset, iBase_FACE);
01517 
01518 #endif
01519   vector<iBase_EntityHandle> entities;
01520   for (unsigned int i = 0; i < nodes.size(); i++)
01521       entities.push_back(nodes[i].gVertexHandle);
01522   iMesh::Error m_err = mk_core()->imesh_instance()->rmvArrTag(&entities[0], entities.size(), m_taghandle);
01523   IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
01524   entities.clear();
01525   for (unsigned int i = 0; i < quads.size(); i++)
01526       entities.push_back(quads[i].gFaceHandle);
01527   m_err = mk_core()->imesh_instance()->rmvArrTag(&entities[0], entities.size(), m_taghandle);
01528   IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
01529   entities.clear();
01530   for (unsigned int i = 0; i < vertices.size(); i++)
01531           entities.push_back(vertices[i].gVertexHandle);
01532   iGeom::Error g_err = mk_core()->igeom_instance()->rmvArrTag(&entities[0], entities.size(), g_taghandle);
01533   IBERRCHK(g_err, "Trouble remove the tag values from an array of entities.");
01534   entities.clear();
01535   for (unsigned int i = 0; i < edges.size(); i++)
01536           entities.push_back(edges[i].gEdgeHandle);
01537   g_err = mk_core()->igeom_instance()->rmvArrTag(&entities[0], entities.size(), g_taghandle);
01538   IBERRCHK(g_err, "Trouble remove the tag values from an array of entities.");
01539 
01540         
01541 }
01542 
01543 //check whether a point is outside the surface or not
01544 bool SubMapping::isOutSideSurf(vector<Vertex> corner, int i, int j)
01545 {//test whether the point (i,j) is inside the polygon or outside
01546    //use the wind number
01547    unsigned int m;
01548    double angle=0;
01549    double p1[2],p2[2];
01550 
01551    for (m = 0; m < corner.size(); m++) {
01552       p1[0] = corner[m].uv[0] - i;
01553       p1[1] = corner[m].uv[1] - j;
01554       p2[0] = corner[(m+1)%corner.size()].uv[0] - i;
01555       p2[1] = corner[(m+1)%corner.size()].uv[1] - j;
01556       angle += Angle2D(p1[0],p1[1],p2[0],p2[1]);
01557    }
01558 
01559    if (fabs(angle) < PI)
01560       return false;
01561    else
01562       return true;
01563 }
01564 
01565 /*
01566    Return the angle between two vectors on a plane
01567    The angle is from vector 1 to vector 2, positive anticlockwise
01568    The result is between -pi -> pi
01569 */
01570 double SubMapping::Angle2D(double x1, double y1, double x2, double y2)
01571 {
01572    double dtheta,theta1,theta2;
01573 
01574    theta1 = atan2(y1,x1);
01575    theta2 = atan2(y2,x2);
01576    dtheta = theta2 - theta1;
01577    while (dtheta > PI)
01578       dtheta -= 2.0*PI;
01579    while (dtheta < -PI)
01580       dtheta += 2.0*PI;
01581 
01582    return dtheta;
01583 }
01584 
01585 }
01586 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines