surfst.cc

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

              for (jt = deleted_triangles.begin();
                   jt != deleted_triangles.end();
                   jt++) {
                  Ref<Triangle> t = *jt;
                  for (j=0; j<2; j++) {
                      Ref<Vertex> v = t->vertex(j);
                      vertices_of_deleted_triangles.insert(v);
                    }
                }

              // find a new point that replaces the deleted vertex
              // (for now use one of the original, since it must lie on the
              // analytic surface)
              Ref<Vertex> replacement_vertex;
              Ref<Edge> short_edge;
              if (hmin==0) {
                  if (l[1] < l[2]) short_edge = tri->edge(1);
                  else short_edge = tri->edge(2);
                }
              else if (hmin==1) {
                  if (l[0] < l[2]) short_edge = tri->edge(0);
                  else short_edge = tri->edge(2);
                }
              else {
                  if (l[0] < l[1]) short_edge = tri->edge(0);
                  else short_edge = tri->edge(1);
                }
              if (short_edge->vertex(0) == vertex) {
                  replacement_vertex = short_edge->vertex(1);
                }
              else {
                  replacement_vertex = short_edge->vertex(0);
                }
              new_vertices.insert(replacement_vertex);
              // for each vertex on the perimeter form a new edge to the
              // replacement vertex (unless the replacement vertex
              // is equal to the perimeter vertex)
              std::map<Ref<Vertex>,Ref<Edge> > new_edge_map;
              for (je = perimeter_edges.begin(); je != perimeter_edges.end();
                   je++) {
                  Ref<Edge> e = *je;
                  for (k = 0; k<2; k++) {
                      Ref<Vertex> v = e->vertex(k);
                      if (v == replacement_vertex) continue;
                      if (new_edge_map.find(v) == new_edge_map.end()) {
                          Ref<Edge> new_e;
                          // if the edge already exists then use the
                          // existing edge
                          for (ke = perimeter_edges.begin();
                               ke != perimeter_edges.end();
                               ke++) {
                              Ref<Edge> tmp = *ke;
                              if ((tmp->vertex(0) == replacement_vertex
                                   &&tmp->vertex(1) == v)
                                  ||(tmp->vertex(1) == replacement_vertex
                                     &&tmp->vertex(0) == v)) {
                                  new_e = tmp;
                                  break;
                                }
                            }
                          if (ke == perimeter_edges.end()) {
                              new_e = newEdge(replacement_vertex,
                                              v);
                            }
                          new_edge_map[v] = new_e;
                          new_edges.insert(new_e);
                        }
                    }
                }
              // for each edge on the perimeter form a new triangle with the
              // replacement vertex (unless the edge contains the replacement
              // vertex)
              for (je = perimeter_edges.begin(); je != perimeter_edges.end();
                   je++) {
                  Ref<Edge> e1 = *je;
                  Ref<Vertex> v0 = e1->vertex(0);
                  Ref<Vertex> v1 = e1->vertex(1);
                  Ref<Edge> e2 = new_edge_map[v0];
                  Ref<Edge> e3 = new_edge_map[v1];
                  if (v0 == replacement_vertex
                      || v1 == replacement_vertex) continue;
                  // Compute the correct orientation of e1 within the new
                  // triangle, by finding the orientation within the old
                  // triangle.
                  int orientation = 0;
                  for (kt = connected_triangles.begin();
                       kt != connected_triangles.end();
                       kt++) {
                      if ((*kt)->contains(e1)) {
                          orientation
                              = (*kt)->orientation(e1);
                          break;
                        }
                    }
                  Ref<Triangle> newtri(newTriangle(e1,e2,e3,orientation));
                  new_triangles.insert(newtri);
                }
            }
        }

#if WRITE_OOGL
      if (_debug) {
          char filename[100];
          static int pass = 0;
          sprintf(filename, "surfst%04d.oogl", pass);
          ExEnv::outn() << scprintf("PASS = %04d\n", pass);
          Ref<Render> render = new OOGLRender(filename);
          Ref<RenderedPolygons> poly = new RenderedPolygons;
          poly->initialize(_vertices.size(), _triangles.size(),
                           RenderedPolygons::Vertex);
          // the number of triangles and edges touching a vertex
          int *n_triangle = new int[_vertices.size()];
          int *n_edge = new int[_vertices.size()];
          memset(n_triangle,0,sizeof(int)*_vertices.size());
          memset(n_edge,0,sizeof(int)*_vertices.size());
          std::set<Ref<Triangle> >::iterator it;
          std::set<Ref<Edge> >::iterator ie;
          std::set<Ref<Vertex> >::iterator iv;
          std::map<Ref<Vertex>, int> vertex_to_index;
          int i = 0;
          for (iv = _vertices.begin(); iv != _vertices.end(); iv++, i++) {
              Ref<Vertex> v = *iv;
              vertex_to_index[v] = i;
              poly->set_vertex(i,
                               v->point()[0],
                               v->point()[1],
                               v->point()[2]);
              if (deleted_vertices.find(v) != deleted_vertices.end()) {
                  poly->set_vertex_rgb(i, 1.0, 0.0, 0.0);
                }
              else {
                  poly->set_vertex_rgb(i, 0.3, 0.3, 0.3);
                }
            }
          i = 0;
          for (it = _triangles.begin(); it != _triangles.end(); it++, i++) {
              Ref<Triangle> t = *it;
              int i0 = vertex_to_index[t->vertex(0)];
              int i1 = vertex_to_index[t->vertex(1)];
              int i2 = vertex_to_index[t->vertex(2)];
              n_triangle[i0]++;
              n_triangle[i1]++;
              n_triangle[i2]++;
              poly->set_face(i,i0,i1,i2);
            }
          for (ie = _edges.begin(); ie != _edges.end(); ie++, i++) {
              Ref<Edge> e = *ie;
              int i0 = vertex_to_index[e->vertex(0)];
              int i1 = vertex_to_index[e->vertex(1)];
              n_edge[i0]++;
              n_edge[i1]++;
            }
          i = 0;
          for (iv = _vertices.begin(); iv != _vertices.end(); iv++, i++) {
              Ref<Vertex> v = *iv;
              if (n_triangle[i] != n_edge[i]) {
                  ExEnv::outn() << "found bad vertex"
                       << " nedge = " << n_edge[i]
                       << " ntriangle = " << n_triangle[i]
                       << endl;
                  if (deleted_vertices.find(v) != deleted_vertices.end()) {
                      poly->set_vertex_rgb(i, 1.0, 1.0, 0.0);
                    }
                  else {
                      poly->set_vertex_rgb(i, 0.0, 1.0, 0.0);
                    }
                }
            }
          render->render(poly.pointer());
          pass++;
          delete[] n_triangle;
          delete[] n_edge;
        }
#endif

      erase_elements_by_value(_triangles,
                              deleted_triangles.begin(),
                              deleted_triangles.end());
      erase_elements_by_value(_edges,
                              deleted_edges.begin(),
                              deleted_edges.end());
      erase_elements_by_value(_vertices,
                              deleted_vertices.begin(),
                              deleted_vertices.end());

      _triangles.insert(new_triangles.begin(), new_triangles.end());
      _edges.insert(new_edges.begin(), new_edges.end());
      _vertices.insert(new_vertices.begin(), new_vertices.end());

      if (_verbose) {
          ExEnv::outn() << "intermediate: ";
          topology_info();
        }

      deleted_edges_length = deleted_edges.size();
    } while(deleted_edges_length != 0);

  // fix the index maps
  _vertex_to_index.clear();
  _edge_to_index.clear();
  _triangle_to_index.clear();

  _index_to_vertex.clear();
  _index_to_edge.clear();
  _index_to_triangle.clear();

  _index_to_vertex.resize(_vertices.size());
  for (i=0, iv = _vertices.begin(); iv != _vertices.end(); i++, iv++) {
      _vertex_to_index[*iv] = i;
      _index_to_vertex[i] = *iv;
    }

  _index_to_edge.resize(_edges.size());
  for (i=0, ie = _edges.begin(); ie != _edges.end(); i++, ie++) {
      _edge_to_index[*ie] = i;
      _index_to_edge[i] = *ie;
    }

  _index_to_triangle.resize(_triangles.size());
  for (i=0, it = _triangles.begin(); it != _triangles.end(); i++, it++) {
      _triangle_to_index[*it] = i;
      _index_to_triangle[i] = *it;
    }

  // restore the int arrays if they were there to begin with
  if (surface_was_completed) {
      complete_int_arrays();
      _completed_surface = 1;
    }

  if (_verbose) {
      ExEnv::outn() << "final: ";
      topology_info();
    }
}


/////////////////////////////////////////////////////////////////////////////

// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: