checktri.cpp

来自「算断裂的」· C++ 代码 · 共 535 行 · 第 1/2 页

          }
          if (notfound) continue;

          if (fdim == sdim) {
            ++numelt1;
            for (int jj = 0; jj < fdim + 1; ++jj) {
              SimpComplex::VertexGlobalIndex vtx = 
                s.node_of_meshface_on_face(fspec, seqno,jj);
              checkoff[vtx] = 1;
            }
          }      
          // Make sure the simplex has the correct orientation.
          
          
          if (checko > 0 && fdim == di) {
            for (int j = 1; j < di + 1; ++j) {
              for (int k = 0; k < di; ++k) {
                simplexmatrix(j - 1, k) = simplexvertex(j,k) - simplexvertex(0,k);
              }
            }
            Real det = simplexmatrix.determinant();
            if ((det <= 0 && checko == 1) || (det >= 0 && checko == 2)) {
              outstream << "?? " << Dumpsimp(s,g,fspec,seqno) 
                <<  " orientation incorrect\n";
              imperfectflag = true;
            }
          }
          
          Real minalt = BIG_REAL;

          // Now find the max side length.
          
          Real maxside = 0.0;
          if (sdim == fdim) {
            for (int i = 0; i < sdim; ++i) {
              for (int j = i + 1; j < sdim + 1; ++j) {
                Real side = 0.0;
                for (int d = 0; d < di; ++d) {
                  Real t = simplexvertex(i, d) - simplexvertex(j,d);
                  side += t*t;
                }
                side = sqrt(side);
                if (side > maxside)
                  maxside = side;
              }
            }
          }
          
          for (int fa = 0; fa < fdim + 1; ++fa) {
            int kk = 0;
            SimpComplex::VertexGlobalIndex lastv;
            for (int jj = 0; jj < fdim + 1; ++jj) {
              if (jj != fa) {
                sortord[kk] = jj;
                key.vert[kk] = s.node_of_meshface_on_face(fspec, seqno, jj);
                ++kk;
              }
              else {
                sortord[fdim] = jj;
                lastv = s.node_of_meshface_on_face(fspec,seqno,jj);
              }
            }
            // Do an insertion sort on the key.     
            int swapcount = 0;
            {
              for (int j = 1; j < fdim; j++) {
                int k1 = j;
                while (k1 > 0 && key.vert[k1] < key.vert[k1-1]) {
                  SimpComplex::VertexGlobalIndex tmp = key.vert[k1];
                  key.vert[k1] = key.vert[k1-1];
                  key.vert[k1-1] = tmp;
                  int tmp2 = sortord[k1];
                  sortord[k1] = sortord[k1 - 1];
                  sortord[k1 - 1] = tmp2;
                  --k1;
                  ++swapcount;
                }
              }
            }
            
            
            // Check the altitude if this simplex is full-dimensional.
            
            // Make a matrix with the data for this simplex facet.
            
            if (fdim == sdim) {
              {
                for (int i = 0; i < sdim; i++) {
                  for (int j = 0; j < di; j++) 
                    simplexmatrix(i,j) = 
                    simplexvertex(sortord[i + 1],j) - simplexvertex(sortord[0],j);
                }
              }
              
              // QR factor the matrix to get the altitude & determinant.
              
              Real det = 1.0; 
              {
                for (int i = 0; i < sdim; i++) {
                  Real beta = Matrix::compute_hh_transform(&hhtrans[i], &simplexmatrix(i,i), 
                    1, di - i);
                  for (int j = i; j < sdim; ++j) {
                    Real ip = 0.0;
                    {
                      for (int k = i; k < di; ++k)
                        ip += simplexmatrix(j,k) * hhtrans[k];
                    }
                    ip *= beta;
                    {
                      for (int k = i; k < di; ++k)
                        simplexmatrix(j,k) -= ip * hhtrans[k];
                    }
                  }
                  det *= simplexmatrix(i,i);
                }
              }    
              
              // Compute the altitude to this facet.
              
              Real alt = fabs(simplexmatrix(sdim - 1,sdim - 1));
              if (alt < minalt)
                minalt = alt;
            }
            
            // hash the facet
            
            Facet_Rec table_entry = facet_occurrence_map[key];     
            if ((swapcount + fa) % 2)
              ++(table_entry.pos_occurrence_count);
            else
              ++(table_entry.neg_occurrence_count);
            int tot_occ = table_entry.pos_occurrence_count + 
              table_entry.neg_occurrence_count;
            if (tot_occ < 3) {
              table_entry.seqno_occ[tot_occ - 1] = seqno;
            }
            facet_occurrence_map[key] = table_entry;
          } // end of facet loop

          if (fdim == sdim) {
            if (maxside > maxsideglobal)
              maxsideglobal = maxside;
            if (minalt < minaltglobal)
              minaltglobal = minalt;        
            if (sdim > 0) {
              Real asp = maxside/minalt;
              if (asp > maxaspglobal) {
                maxaspglobal = asp;
                worsttriseq = seqno;
                worsttrifspec = fspec;
              }
            }
          }
          // Done with this element
        }
      }
      // Now go over the facets to make sure that each facet occurs the proper
      // number of times.

      // First, loop over the facets that belong to lower-dimensional
      // topological entities.

      for (Brep::Face_Spec_Loop_Over_Face_Subfaces subfspec(g, fspec);
      subfspec.notdone();
      ++subfspec) {
        if (subfspec.fdim() != fdim - 1) continue;
        int nmf2 = s.num_meshface_on_face(subfspec);
        Facet_Key key;
        for (int seqno2 = 0; seqno2 < nmf2; ++seqno2) {
          for (int jj = 0; jj < fdim; ++jj)
            key.vert[jj] = s.node_of_meshface_on_face(subfspec, seqno2, jj);
          // Insertion sort on key.          
          {
            for (int j = 1; j < fdim; ++j) {
              int k1 = j;
              while (k1 > 0 && key.vert[k1] < key.vert[k1-1]) {
                SimpComplex::VertexGlobalIndex tmp = key.vert[k1];
                key.vert[k1] = key.vert[k1-1];
                key.vert[k1-1] = tmp;
                --k1;
              }
            }
          }
          Facet_Rec table_entry = facet_occurrence_map[key];
          ++table_entry.occurs_on_bdry;
          facet_occurrence_map[key] = table_entry;
        }
      }

      // Next pass: make sure counts for boundaries add up correctly.

      for (Facet_Map_Type::const_iterator it = facet_occurrence_map.begin();
      it != facet_occurrence_map.end();
      ++it) {
        Facet_Key key = it -> first;
        Facet_Rec table_entry = it -> second;
        int pos1 = table_entry.pos_occurrence_count;
        int neg1 = table_entry.neg_occurrence_count;
        int tot1 = pos1 + neg1;
        int bdry1 = table_entry.occurs_on_bdry;
        const char* msg = "";
        if (bdry1 > 0 && tot1 == 0) {
          msg = "mesh entity occur on topological boundary but not as boundary of mesh";
        }
        else if (bdry1 > 0 && ((tot1 + bdry1) % 2))  {
          msg = "mesh entity occurs on topological boundary with wrong parity";
        }
        else if (bdry1 == 0 && tot1 % 2) {
          msg = "internal mesh entity occurs an odd number of times (hole in mesh)";
        }
        if (fdim >= di - 1 && (pos1 > 1 || neg1 > 1)) 
          msg = "mesh entity contained twice on the same side";
        if (strlen(msg) > 0) {
          outstream << "?? Problem with mesh entity bounded by nodes ";
          for (int j = 0; j < fdim; ++j)
            outstream << key.vert[j] << " ";
          outstream << " " << msg 
           << "\n This mesh entity is contained by the following higher dimensional\n"
           << " mesh entities:\n";
          int nume1 = (tot1 > 2)? 2 : tot1;
          for (int j2 = 0; j2 < nume1; ++j2)
            outstream << Dumpsimp(s,g,fspec,table_entry.seqno_occ[j2]) << '\n';
          if (tot1 > 2)
            outstream << " and others.\n";
          imperfectflag = true;
        }
      }
    }
  }

  // Make sure every vertex is checked off.

  {
    for (SimpComplex::VertexOrdinalIndex vnumo = 0; 
    vnumo < numvtx;
    ++vnumo) {
      SimpComplex::VertexGlobalIndex vnum = s.ordinal_to_global(vnumo);
      if (checkoff[vnum] == 0) {
        outstream << "?? Vertex " << vnum << " not used in any simplex\n";
        imperfectflag = true;
      }
    }
  }

  int numelt = s.num_elements();
  if (numelt1 != numelt) {
    outstream << "?? Mismatch: stated number of elements = " << numelt
      << " versus actual number = " << numelt1 << '\n';
    imperfectflag = true;
  }
      

  if (worsttriseq >= 0) {
    outstream << "Maximum aspect ratio =        " << maxaspglobal 
      << " achieved in\n" << Dumpsimp(s,g,worsttrifspec,worsttriseq);
    outstream << "Maximum global side length =  " << maxsideglobal << "\n";
    outstream << "Minimum global altitude =     " << minaltglobal << "\n";
  }
  outstream << "Number of nodes = " << s.num_nodes() 
	    << " number of elements = " << s.num_elements() << '\n';
  if (imperfectflag) {
   outstream << "!! One or more flaws were detected in the mesh\n";
   maxaspglobal = -1.0;
  }
  return maxaspglobal;
}