surf.cc

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

//
// surf.cc
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <[e-mail]>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//

#ifdef __GNUC__
#pragma implementation
#endif

#include <util/misc/formio.h>
#include <util/keyval/keyval.h>
#include <math/scmat/matrix.h>
#include <math/scmat/vector3.h>
#include <math/isosurf/surf.h>
#include <math/isosurf/isosurf.h>
#include <util/render/polygons.h>

using namespace std;
using namespace sc;

#ifndef WRITE_OOGL
#define WRITE_OOGL 1
#endif

#if WRITE_OOGL
#include <util/render/oogl.h>
#endif

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

TriangulatedSurface::TriangulatedSurface():
  _verbose(0),
  _debug(0),
  _triangle_vertex(0),
  _triangle_edge(0),
  _edge_vertex(0),
  _integrator(new GaussTriangleIntegrator(1))
{
  clear();
}

TriangulatedSurface::TriangulatedSurface(const Ref<KeyVal>& keyval):
  _triangle_vertex(0),
  _triangle_edge(0),
  _edge_vertex(0)
{
  _verbose = keyval->booleanvalue("verbose");
  _debug = keyval->booleanvalue("debug");
  Ref<TriangleIntegrator> triint;
  triint << keyval->describedclassvalue("integrator");
  if (triint.null()) {
      triint = new GaussTriangleIntegrator(1);
    }
  set_integrator(triint);
  triint << keyval->describedclassvalue("fast_integrator");
  set_fast_integrator(triint);
  triint << keyval->describedclassvalue("accurate_integrator");
  set_accurate_integrator(triint);
  clear();
}

TriangulatedSurface::~TriangulatedSurface()
{
  clear();
}

void
TriangulatedSurface::topology_info(ostream&o)
{
  topology_info(nvertex(), nedge(), ntriangle(), o);
}

void
TriangulatedSurface::topology_info(int v, int e, int t, ostream&o)
{
  // Given v vertices i expect 2*v - 4*n_surface triangles
  // and 3*v - 6*n_surface edges
  o << indent
    << scprintf("n_vertex = %d, n_edge = %d, n_triangle = %d:",
                v, e, t)
    << endl;
  int nsurf_e = ((3*v - e)%6 == 0)? (3*v - e)/6 : -1;
  int nsurf_t = ((2*v - t)%4 == 0)? (2*v - t)/4 : -1;
  if ((nsurf_e!=-1) && (nsurf_e == nsurf_t)) {
      o << indent
        << scprintf("  this is consistent with n_closed_surface - n_hole = %d",
                    nsurf_e)
        << endl;
    }
  else {
      o << indent
        << scprintf("  this implies that some surfaces are not closed")
        << endl;
    }
}

void
TriangulatedSurface::set_integrator(const Ref<TriangleIntegrator>& i)
{
  _integrator = i;
}

void
TriangulatedSurface::set_fast_integrator(const Ref<TriangleIntegrator>& i)
{
  _fast_integrator = i;
}

void
TriangulatedSurface::set_accurate_integrator(const Ref<TriangleIntegrator>& i)
{
  _accurate_integrator = i;
}

Ref<TriangleIntegrator>
TriangulatedSurface::integrator(int)
{
  // currently the argument, the integer index of the triangle, is ignored
  return _integrator;
}

Ref<TriangleIntegrator>
TriangulatedSurface::fast_integrator(int)
{
  // currently the argument, the integer index of the triangle, is ignored
  return _fast_integrator.null()?_integrator:_fast_integrator;
}

Ref<TriangleIntegrator>
TriangulatedSurface::accurate_integrator(int)
{
  // currently the argument, the integer index of the triangle, is ignored
  return _accurate_integrator.null()?_integrator:_accurate_integrator;
}

void
TriangulatedSurface::clear_int_arrays()
{
  if (_triangle_vertex) {
      for (int i=0; i<_triangles.size(); i++) {
          delete[] _triangle_vertex[i];
        }
      delete[] _triangle_vertex;
    }
  _triangle_vertex = 0;

  if (_triangle_edge) {
      for (int i=0; i<_triangles.size(); i++) {
          delete[] _triangle_edge[i];
        }
      delete[] _triangle_edge;
    }
  _triangle_edge = 0;

  if (_edge_vertex) {
      for (int i=0; i<_edges.size(); i++) {
          delete[] _edge_vertex[i];
        }
      delete[] _edge_vertex;
    }
  _edge_vertex = 0;

  _completed_surface = 0;
}

void
TriangulatedSurface::clear()
{
  _completed_surface = 0;

  clear_int_arrays();

  _have_values = 0;
  _values.clear();

  _vertices.clear();
  _edges.clear();
  _triangles.clear();

  _tmp_edges.clear();
}

void
TriangulatedSurface::complete_surface()
{
  complete_ref_arrays();
  complete_int_arrays();

  _completed_surface = 1;
}

void
TriangulatedSurface::complete_ref_arrays()
{
  _tmp_edges.clear();
  _index_to_edge.clear();
  _edge_to_index.clear();

  int i;
  int ntri = ntriangle();
  _edges.clear();
  for (i=0; i<ntri; i++) {
      Ref<Triangle> tri = triangle(i);
      add_edge(tri->edge(0));
      add_edge(tri->edge(1));
      add_edge(tri->edge(2));
    }
  int ne = nedge();
  _vertices.clear();
  _index_to_vertex.clear();
  _vertex_to_index.clear();
  for (i=0; i<ne; i++) {
      Ref<Edge> e = edge(i);
      add_vertex(e->vertex(0));
      add_vertex(e->vertex(1));
    }
}

void
TriangulatedSurface::complete_int_arrays()
{
  clear_int_arrays();
  
  int i;
  int ntri = ntriangle();
  int ne = nedge();

  // construct the array that converts the triangle number and vertex
  // number within the triangle to the overall vertex number
  _triangle_vertex = new int*[ntri];
  for (i=0; i<ntri; i++) {
      _triangle_vertex[i] = new int[3];
      for (int j=0; j<3; j++) {
          Ref<Vertex> v = triangle(i)->vertex(j);
          _triangle_vertex[i][j] = _vertex_to_index[v];
        }
    }

  // construct the array that converts the triangle number and edge number
  // within the triangle to the overall edge number
  _triangle_edge = new int*[ntri];
  for (i=0; i<ntri; i++) {
      _triangle_edge[i] = new int[3];
      for (int j=0; j<3; j++) {
          Ref<Edge> e = triangle(i)->edge(j);
          _triangle_edge[i][j] = _edge_to_index[e];
        }
    }

  // construct the array that converts the edge number and vertex number
  // within the edge to the overall vertex number
  _edge_vertex = new int*[ne];
  for (i=0; i<ne; i++) {
      _edge_vertex[i] = new int[2];
      for (int j=0; j<2; j++) {
          Ref<Vertex> v = edge(i)->vertex(j);
          _edge_vertex[i][j] = _vertex_to_index[v];
        }
    }
}

void
TriangulatedSurface::compute_values(Ref<Volume>&vol)
{
  int n = _vertices.size();
  _values.resize(n);

  for (int i=0; i<n; i++) {
      vol->set_x(vertex(i)->point());
      _values[i] = vol->value();
    }
  _have_values = 1;
}

double
TriangulatedSurface::flat_area()
{
  double result = 0.0;
  for (std::set<Ref<Triangle> >::iterator i=_triangles.begin();
       i!=_triangles.end(); i++) {
      result += (*i)->flat_area();
    }
  return result;
}

double
TriangulatedSurface::flat_volume()
{
  double result = 0.0;
  for (int i=0; i<_triangles.size(); i++) {

      // get the vertices of the triangle
      SCVector3 A(vertex(triangle_vertex(i,0))->point());
      SCVector3 B(vertex(triangle_vertex(i,1))->point());
      SCVector3 C(vertex(triangle_vertex(i,2))->point());

      // project the vertices onto the xy plane
      SCVector3 Axy(A); Axy[2] = 0.0;
      SCVector3 Bxy(B); Bxy[2] = 0.0;
      SCVector3 Cxy(C); Cxy[2] = 0.0;

      // construct the legs of the triangle in the xy plane
      SCVector3 BAxy = Bxy - Axy;
      SCVector3 CAxy = Cxy - Axy;

      // find the lengths of the legs of the triangle in the xy plane
      double baxy = sqrt(BAxy.dot(BAxy));
      double caxy = sqrt(CAxy.dot(CAxy));

      // if one of the legs is of length zero, then there is
      // no contribution from this triangle
      if (baxy < 1.e-16 || caxy < 1.e-16) continue;

      // find the sine of the angle between the legs of the triangle
      // in the xy plane
      double costheta = BAxy.dot(CAxy)/(baxy*caxy);
      double sintheta = sqrt(1.0 - costheta*costheta);

      // the area of the triangle in the xy plane
      double areaxy = 0.5 * baxy * caxy * sintheta;

      // the height of the three corners of the triangle
      // (relative to the z plane)
      double hA = A[2];
      double hB = B[2];
      double hC = C[2];

      // the volume of the space under the triangle
      double volume = areaxy * (hA + (hB + hC - 2.0*hA)/3.0);

      // the orientation of the triangle along the projection axis (z)
      SCVector3 BA(B-A);
      SCVector3 CA(C-A);
      double z_orientation = BA.cross(CA)[2];

      if (z_orientation > 0.0) {
          result += volume;
        }
      else {
          result -= volume;
        }

    }

  // If the volume is negative, then the surface gradients were
  // opposite in sign to the direction assumed.  Flip the sign
  // to fix.
  return fabs(result);
}

double
TriangulatedSurface::area()
{
  double area = 0.0;
  TriangulatedSurfaceIntegrator triint(this);
  for (triint = 0; triint.update(); triint++) {
      area += triint.w();
    }
  return area;
}

double
TriangulatedSurface::volume()
{
  double volume = 0.0;
  TriangulatedSurfaceIntegrator triint(this);
  for (triint = 0; triint.update(); triint++) {
      volume += triint.weight()*triint.dA()[2]*triint.current()->point()[2];
    }
  return volume;
}

void
TriangulatedSurface::add_vertex(const Ref<Vertex>&t)
{
  int i = _vertices.size();
  _vertices.insert(t);
  if (i != _vertices.size()) {
      _index_to_vertex.push_back(t);
      _vertex_to_index[t] = i;
      if (_index_to_vertex.size() != _vertex_to_index.size()) {
          ExEnv::errn() << "TriangulatedSurface::add_vertex: length mismatch" << endl;
          abort();
        }
    }
}

void
TriangulatedSurface::add_edge(const Ref<Edge>&t)
{
  int i = _edges.size();
  _edges.insert(t);
  if (i != _edges.size()) {
      _index_to_edge.push_back(t);
      _edge_to_index[t] = i;
      if (_index_to_edge.size() != _edge_to_index.size()) {
          ExEnv::errn() << "TriangulatedSurface::add_edge: length mismatch" << endl;
          abort();
        }
    }
}

void
TriangulatedSurface::add_triangle(const Ref<Triangle>&t)
{
  if (_completed_surface) clear();
  int i = _triangles.size();
  _triangles.insert(t);
  if (i != _triangles.size()) {
      _index_to_triangle.push_back(t);
      _triangle_to_index[t] = i;
      if (_index_to_triangle.size() != _triangle_to_index.size()) {
          ExEnv::errn() << "TriangulatedSurface::add_triangle: length mismatch" << endl;
          abort();
        }
    }
}

void
TriangulatedSurface::add_triangle(const Ref<Vertex>& v1,
                                  const Ref<Vertex>& v2,
                                  const Ref<Vertex>& v3)
{
  // Find this triangle's edges if they have already be created
  // for some other triangle.
  Ref<Edge> e0, e1, e2;

  const std::set<Ref<Edge> > &v1edges = _tmp_edges[v1];

  const std::set<Ref<Edge> > &v2edges = _tmp_edges[v2];

  std::set<Ref<Edge> >::const_iterator ix;
  for (ix = v1edges.begin(); ix != v1edges.end(); ix++) {
      const Ref<Edge>& e = *ix;
      if (e->vertex(0) == v2 || e->vertex(1) == v2) {
          e0 = e;
        }
      else if (e->vertex(0) == v3 || e->vertex(1) == v3) {
          e2 = e;
        }
    }
  for (ix = v2edges.begin(); ix != v2edges.end(); ix++) {
      const Ref<Edge>& e = *ix;
      if (e->vertex(0) == v3 || e->vertex(1) == v3) {
          e1 = e;
        }
    }

  if (e0.null()) {
      e0 = newEdge(v1,v2);
      _tmp_edges[v1].insert(e0);
      _tmp_edges[v2].insert(e0);
    }
  if (e1.null()) {
      e1 = newEdge(v2,v3);
      _tmp_edges[v2].insert(e1);
      _tmp_edges[v3].insert(e1);
    }
  if (e2.null()) {
      e2 = newEdge(v3,v1);