surf.cc

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

      _tmp_edges[v3].insert(e2);
      _tmp_edges[v1].insert(e2);
    }
  
  int orientation;
  if (e0->vertex(0) == v1) {
      orientation = 0;
    }
  else {
      orientation = 1;
    }
  
  add_triangle(newTriangle(e0,e1,e2,orientation));
}

// If a user isn't keeping track of edges while add_triangle is being
// used to build the surface, then this can be called to see if an edge
// already exists (at a great performance cost).
Ref<Edge>
TriangulatedSurface::find_edge(const Ref<Vertex>& v1, const Ref<Vertex>& v2)
{
  std::set<Ref<Triangle> >::iterator i;

  for (i=_triangles.begin(); i!=_triangles.end(); i++) {
      Ref<Triangle> t = *i;
      Ref<Edge> e1 = t->edge(0);
      Ref<Edge> e2 = t->edge(1);
      Ref<Edge> e3 = t->edge(2);
      if (e1->vertex(0) == v1 && e1->vertex(1) == v2) return e1;
      if (e1->vertex(1) == v1 && e1->vertex(0) == v2) return e1;
      if (e2->vertex(0) == v1 && e2->vertex(1) == v2) return e2;
      if (e2->vertex(1) == v1 && e2->vertex(0) == v2) return e2;
      if (e3->vertex(0) == v1 && e3->vertex(1) == v2) return e3;
      if (e3->vertex(1) == v1 && e3->vertex(0) == v2) return e3;
    }

  return 0;
}

void
TriangulatedSurface::print(ostream&o) const
{
  o << indent << "TriangulatedSurface:" << endl;
  int i;

  int np = nvertex();
  o << indent << scprintf(" %3d Vertices:",np) << endl;
  for (i=0; i<np; i++) {
      Ref<Vertex> p = vertex(i);
      o << indent << scprintf("  %3d:",i);
      for (int j=0; j<3; j++) {
          o << scprintf(" % 15.10f", p->point()[j]);
        }
      o << endl;
    }

  int ne = nedge();
  o << indent << scprintf(" %3d Edges:",ne) << endl;
  for (i=0; i<ne; i++) {
      Ref<Edge> e = edge(i);
      Ref<Vertex> v0 = e->vertex(0);
      Ref<Vertex> v1 = e->vertex(1);
      std::map<Ref<Vertex>,int>::const_iterator v0i=_vertex_to_index.find(v0);
      std::map<Ref<Vertex>,int>::const_iterator v1i=_vertex_to_index.find(v1);
      int v0int = v0i==_vertex_to_index.end()? -1: v0i->second;
      int v1int = v1i==_vertex_to_index.end()? -1: v1i->second;
      o << indent
        << scprintf("  %3d: %3d %3d",i, v0int, v1int)
        << endl;
    }

  int nt = ntriangle();
  o << indent << scprintf(" %3d Triangles:",nt) << endl;
  for (i=0; i<nt; i++) {
      Ref<Triangle> tri = triangle(i);
      Ref<Edge> e0 = tri->edge(0);
      Ref<Edge> e1 = tri->edge(1);
      Ref<Edge> e2 = tri->edge(2);
      std::map<Ref<Edge>,int>::const_iterator e0i = _edge_to_index.find(e0);
      std::map<Ref<Edge>,int>::const_iterator e1i = _edge_to_index.find(e1);
      std::map<Ref<Edge>,int>::const_iterator e2i = _edge_to_index.find(e2);
      int e0int = e0i==_edge_to_index.end()? -1: e0i->second;
      int e1int = e1i==_edge_to_index.end()? -1: e1i->second;
      int e2int = e2i==_edge_to_index.end()? -1: e2i->second;
      o << indent
        << scprintf("  %3d: %3d %3d %3d",i, e0int, e1int, e2int)
        << endl;
    }
}

void
TriangulatedSurface::print_vertices_and_triangles(ostream&o) const
{
  o << indent << "TriangulatedSurface:" << endl;
  int i;

  int np = nvertex();
  o << indent << scprintf(" %3d Vertices:",np) << endl;
  for (i=0; i<np; i++) {
      Ref<Vertex> p = vertex(i);
      o << indent << scprintf("  %3d:",i);
      for (int j=0; j<3; j++) {
          o << scprintf(" % 15.10f", p->point()[j]);
        }
      o << endl;
    }

  int nt = ntriangle();
  o << indent << scprintf(" %3d Triangles:",nt) << endl;
  for (i=0; i<nt; i++) {
      Ref<Triangle> tri = triangle(i);
      o << indent
        << scprintf("  %3d: %3d %3d %3d",i,
                    _triangle_vertex[i][0],
                    _triangle_vertex[i][1],
                    _triangle_vertex[i][2])
        << endl;
    }
}

void
TriangulatedSurface::render(const Ref<Render> &render)
{
  Ref<RenderedPolygons> poly = new RenderedPolygons;
  poly->initialize(_vertices.size(), _triangles.size(),
                   RenderedPolygons::Vertex);
  std::set<Ref<Vertex> >::iterator iv;
  std::set<Ref<Triangle> >::iterator it;
  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]);
      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;
      poly->set_face(i,
                     vertex_to_index[t->vertex(0)],
                     vertex_to_index[t->vertex(1)],
                     vertex_to_index[t->vertex(2)]);
    }
  render->render(poly.pointer());
}

void
TriangulatedSurface::print_geomview_format(ostream&o) const
{
  o << "OFF" << endl;

  o << nvertex() << " " << ntriangle() << " " << nedge() << endl;
  int i;

  int np = nvertex();
  for (i=0; i<np; i++) {
      Ref<Vertex> p = vertex(i);
      for (int j=0; j<3; j++) {
          o << scprintf(" % 15.10f", p->point()[j]);
        }
      o << endl;
    }

  int nt = ntriangle();
  for (i=0; i<nt; i++) {
      Ref<Triangle> tri = triangle(i);
      o << scprintf(" 3 %3d %3d %3d",
                    _triangle_vertex[i][0],
                    _triangle_vertex[i][1],
                    _triangle_vertex[i][2])
          << endl;
    }
}

void
TriangulatedSurface::recompute_index_maps()
{
  int i;
  std::set<Ref<Vertex> >::iterator iv;
  std::set<Ref<Edge> >::iterator ie;
  std::set<Ref<Triangle> >::iterator it;

  // 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;
    }

}

Edge*
TriangulatedSurface::newEdge(const Ref<Vertex>& v0, const Ref<Vertex>& v1) const
{
  return new Edge(v0,v1);
}

Triangle*
TriangulatedSurface::newTriangle(const Ref<Edge>& e0,
                                 const Ref<Edge>& e1,
                                 const Ref<Edge>& e2,
                                 int orientation) const
{
  return new Triangle(e0,e1,e2,orientation);
}

//////////////////////////////////////////////////////////////////////
// TriangulatedSurfaceIntegrator

TriangulatedSurfaceIntegrator::
  TriangulatedSurfaceIntegrator(const Ref<TriangulatedSurface>&ts)
{
  set_surface(ts);
  use_default_integrator();

  _itri = 0;
  _irs = 0;
}

TriangulatedSurfaceIntegrator::
  TriangulatedSurfaceIntegrator()
{
  use_default_integrator();
}

void
TriangulatedSurfaceIntegrator::
  operator =(const TriangulatedSurfaceIntegrator&i)
{
  set_surface(i._ts);

  _integrator = i._integrator;
  _itri = i._itri;
  _irs = i._irs;
}

TriangulatedSurfaceIntegrator::
  ~TriangulatedSurfaceIntegrator()
{
}

void
TriangulatedSurfaceIntegrator::set_surface(const Ref<TriangulatedSurface>&s)
{
  _ts = s;
  _current = new Vertex();
}

int
TriangulatedSurfaceIntegrator::
  vertex_number(int i)
{
  return _ts->triangle_vertex(_itri,i);
}

Ref<Vertex>
TriangulatedSurfaceIntegrator::
  current()
{
  return _current;
}

int
TriangulatedSurfaceIntegrator::n()
{
  int result = 0;
  int ntri = _ts->ntriangle();
  for (int i=0; i<ntri; i++) {
      result += (_ts.pointer()->*_integrator)(i)->n();
    }
  return result;
}

int
TriangulatedSurfaceIntegrator::update()
{
  if (_itri < 0 || _itri >= _ts->ntriangle()) return 0;

  TriangleIntegrator* i = (_ts.pointer()->*_integrator)(_itri).pointer();
  _s = i->s(_irs);
  _r = i->r(_irs);
  _weight = i->w(_irs);
  Ref<Triangle> t = _ts->triangle(_itri);
  Ref<TriInterpCoef> coef = i->coef(t->order(),_irs);
  t->interpolate(coef, _r, _s, _current, _dA);
  _surface_element = _dA.norm();
  static double cum;
  if (_irs == 0) cum = 0.0;
  cum += _surface_element * _weight;
  //ExEnv::outn() << scprintf("%2d dA = %12.8f, w = %12.8f, Sum wdA = %12.8f",
  //                 _irs, _surface_element, _weight, cum)
  //     << endl;

  return (int) 1;
}

void
TriangulatedSurfaceIntegrator::
  operator ++()
{
  int n = (_ts.pointer()->*_integrator)(_itri)->n();
  if (_irs == n-1) {
      _irs = 0;
      if (_grp.null()) _itri++;
      else _itri += _grp->n();
    }
  else {
      _irs++;
    }
}

void
TriangulatedSurfaceIntegrator::distribute(const Ref<MessageGrp> &grp)
{
  _grp = grp;
}

int
TriangulatedSurfaceIntegrator::
  operator = (int i)
{
  _itri = i;
  _irs = 0;
  return i;
}

void
TriangulatedSurfaceIntegrator::use_default_integrator()
{
  _integrator = &TriangulatedSurface::integrator;
}

void
TriangulatedSurfaceIntegrator::use_fast_integrator()
{
  _integrator = &TriangulatedSurface::fast_integrator;
}

void
TriangulatedSurfaceIntegrator::use_accurate_integrator()
{
  _integrator = &TriangulatedSurface::accurate_integrator;
}

/////////////////////////////////////////////////////////////////////////
// TriangulatedImplicitSurface
static ClassDesc TriangulatedImplicitSurface_cd(
  typeid(TriangulatedImplicitSurface),"TriangulatedImplicitSurface",1,"public TriangulatedSurface",
  0, create<TriangulatedImplicitSurface>, 0);

TriangulatedImplicitSurface::
TriangulatedImplicitSurface(const Ref<KeyVal>&keyval):
  TriangulatedSurface(keyval)
{
  inited_ = 0;

  vol_ << keyval->describedclassvalue("volume");
  if (keyval->error() != KeyVal::OK) {
      ExEnv::errn() << "TriangulatedImplicitSurface(const Ref<KeyVal>&keyval): "
           << "requires \"volume\"" << endl;
      abort();
    }

  isovalue_ = keyval->doublevalue("value");
  if (keyval->error() != KeyVal::OK) isovalue_ = 0.0;

  fix_orientation_ = keyval->booleanvalue("fix_orientation");
  if (keyval->error() != KeyVal::OK) fix_orientation_ = 1;

  remove_short_edges_ = keyval->booleanvalue("remove_short_edges");
  if (keyval->error() != KeyVal::OK) remove_short_edges_ = 1;

  remove_slender_triangles_ = keyval->booleanvalue("remove_slender_triangles");
  if (keyval->error() != KeyVal::OK) remove_slender_triangles_ = 0;

  remove_small_triangles_ = keyval->booleanvalue("remove_small_triangles");
  if (keyval->error() != KeyVal::OK) remove_small_triangles_ = 0;

  short_edge_factor_ = keyval->doublevalue("short_edge_factor");
  if (keyval->error() != KeyVal::OK) short_edge_factor_ = 0.4;

  slender_triangle_factor_ = keyval->doublevalue("slender_triangle_factor");
  if (keyval->error() != KeyVal::OK) slender_triangle_factor_ = 0.2;

  small_triangle_factor_ = keyval->doublevalue("small_triangle_factor");
  if (keyval->error() != KeyVal::OK) small_triangle_factor_ = 0.2;

  resolution_ = keyval->doublevalue("resolution");
  if (keyval->error() != KeyVal::OK) resolution_ = 1.0;

  order_ = keyval->intvalue("order");
  if (keyval->error() != KeyVal::OK) order_ = 1;

  int initialize = keyval->booleanvalue("initialize");
  if (initialize) init();
}

void
TriangulatedImplicitSurface::init()
{
  if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: init start" << endl;
  ImplicitSurfacePolygonizer isogen(vol_);
  isogen.set_resolution(resolution_);

  if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: isosurface" << endl;
  isogen.isosurface(isovalue_,*this);
#if WRITE_OOGL
  if (_debug) {
      render(new OOGLRender("surfiso.oogl"));
    }
#endif
  if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: orientation" << endl;
  if (fix_orientation_) fix_orientation();
#if WRITE_OOGL
  if (_debug) {
      render(new OOGLRender("surffix.oogl"));
    }
#endif
  if (remove_short_edges_) {
      if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: short edges" << endl;
      remove_short_edges(short_edge_factor_*resolution_,vol_,isovalue_);
      if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: orientation" << endl;
      if (fix_orientation_) fix_orientation();
    }
  if (remove_slender_triangles_ || remove_small_triangles_) {
      if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: slender" << endl;
      double height_cutoff = slender_triangle_factor_ * resolution_;
      double area_cutoff = small_triangle_factor_*resolution_*resolution_*0.5;
      remove_slender_triangles(remove_slender_triangles_, height_cutoff,
                               remove_small_triangles_, area_cutoff,
                               vol_,isovalue_);
      if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: orientation" << endl;
      if (fix_orientation_) fix_orientation();
    }

  // see if a higher order approximation to the surface is required
  if (order_ > 1) {
      if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: order" << endl;
      int i;
      for (i=0; i<nedge(); i++) {
          edge(i)->set_order(order_, vol_, isovalue_);
        }
      for (i=0; i<ntriangle(); i++) {
          triangle(i)->set_order(order_, vol_, isovalue_);
        }
    }
  inited_ = 1;
  if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: init done" << endl;
}

TriangulatedImplicitSurface::~TriangulatedImplicitSurface()
{
}

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

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