triangle.cc

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

//
// triangle.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/isosurf/triangle.h>
#include <math/scmat/vector3.h>

using namespace std;
using namespace sc;

/////////////////////////////////////////////////////////////////////////
// Triangle
///////////////////////////////////////////////////////////////////////////
// Here is the layout of the vertices and edges in the triangles:         |
//                                                                        |
//                       vertex(1) (r=0, s=1)                             |
//                          +                                             |
//                         / \  _edges[1].vertex(_orientation1)           |
//                        /   \                                           |
//                       /     \                                          |
//                      /       \                                         |
//                     /         \                                        |
//                    /           \                                       |
//         _edges[0] /             \  _edges[1]                           |
//          (r = 0) /               \   (1-r-s = 0)                       |
//                 /                 \                                    |
//                /                   \                                   |
//               /                     \                                  |
//              /                       \ _edges[1].vertex(!_orientation1)|
//             /                         \                                |
//   vertex(0)+---------------------------+                               |
// (r=0, s=0)        _edges[2] (s = 0)      vertex(2) (r=1, s=0)          |
//                                                                        |
//  Zienkiewicz and Taylor, "The Finite Element Method", 4th Ed, Vol 1,   |
//  use                                                                   |
//      L1 = 1 - r - s                                                    |
//      L2 = r,                                                           |
//      L3 = s.                                                           |
//  I also use these below.                                               |
///////////////////////////////////////////////////////////////////////////

Triangle::Triangle(const Ref<Edge>& v1, const Ref<Edge>& v2, const Ref<Edge>& v3,
                   unsigned int orientation0)
{
  _orientation0 = orientation0;

  _edges[0] = v1;
  _edges[1] = v2;
  _edges[2] = v3;

  // edge 0 corresponds to r = 0
  // edge 1 corresponds to (1-r-s) = 0
  // edge 2 corresponds to s = 0
  // edge 0 vertex _orientation0 is (r=0,s=0)
  // edge 1 vertex _orientation1 is (r=0,s=1)
  // edge 2 vertex _orientation2 is (r=1,s=0)
  // edge 0 vertex (1-_orientation0) is edge 1 vertex _orientation1
  // edge 1 vertex (1-_orientation1) is edge 2 vertex _orientation2
  // edge 2 vertex (1-_orientation2) is edge 0 vertex _orientation0

  Ref<Edge> *e = _edges;

  // swap edges 1 and 2 if necessary
  if (e[0]->vertex(1-_orientation0) != e[1]->vertex(0)) {
      if (e[0]->vertex(1-_orientation0) != e[1]->vertex(1)) {
          e[1] = v3;
          e[2] = v2;
        }
    }

  // compute the orientation of _edge[1]
  if (e[0]->vertex(1-_orientation0) == e[1]->vertex(0)) {
      _orientation1 = 0;
    }
  else {
      _orientation1 = 1;
    }

  // compute the orientation of _edge[2]
  if (e[1]->vertex(1-_orientation1) == e[2]->vertex(0)) {
      _orientation2 = 0;
    }
  else {
      _orientation2 = 1;
    }

  // make sure that the triangle is really a triangle
  if ( e[0]->vertex(1-_orientation0) != e[1]->vertex(_orientation1)
       || e[1]->vertex(1-_orientation1) != e[2]->vertex(_orientation2)
       || e[2]->vertex(1-_orientation2) != e[0]->vertex(_orientation0))
    {
      ExEnv::errn() << "Triangle: given edges that don't form a triangle" << endl;
      abort();
    }

  _order = 1;
  _vertices = new Ref<Vertex>[3];

  _vertices[TriInterpCoef::ijk_to_index(_order, 0, 0)] = vertex(0);
  _vertices[TriInterpCoef::ijk_to_index(0, 0, _order)] = vertex(1);
  _vertices[TriInterpCoef::ijk_to_index(0, _order, 0)] = vertex(2);
}

Triangle::~Triangle()
{
  if (_vertices) delete[] _vertices;
}

Ref<Vertex>
Triangle::vertex(int i)
{
  return _edges[i]->vertex(orientation(i));
}

unsigned int
Triangle::orientation(const Ref<Edge>& e) const
{
  if (e == _edges[0]) return orientation(0);
  if (e == _edges[1]) return orientation(1);
  return orientation(2);
}

int
Triangle::contains(const Ref<Edge>& e) const
{
  if (_edges[0] == e) return 1;
  if (_edges[1] == e) return 1;
  if (_edges[2] == e) return 1;
  return 0;
}

double Triangle::flat_area()
{
  double a = _edges[0]->straight_length();
  double b = _edges[1]->straight_length();
  double c = _edges[2]->straight_length();
  double a2 = a*a;
  double b2 = b*b;
  double c2 = c*c;
  return 0.25 * sqrt( 2.0 * (c2*b2 + c2*a2 + a2*b2)
                      - c2*c2 - b2*b2 - a2*a2);
}

void Triangle::add_vertices(std::set<Ref<Vertex> >&set)
{
  for (int i=0; i<3; i++) set.insert(_edges[i]->vertex(orientation(i)));
}

void Triangle::add_edges(std::set<Ref<Edge> >&set)
{
  for (int i=0; i<3; i++) set.insert(_edges[i]);
}

void
Triangle::interpolate(double r,double s,const Ref<Vertex>&result,SCVector3&dA)
{
  TriInterpCoefKey key(_order, r, s);
  Ref<TriInterpCoef> coef = new TriInterpCoef(key);
  interpolate(coef, r, s, result, dA);
}

void
Triangle::interpolate(const Ref<TriInterpCoef>& coef,
                      double r, double s,
                      const Ref<Vertex>&result, SCVector3&dA)
{
  unsigned int i, j, k;

  //double L1 = 1 - r - s;
  //double L2 = r;
  //double L3 = s;

  SCVector3 tmp(0.0);
  SCVector3 x_s(0.0);
  SCVector3 x_r(0.0);
  for (i=0; i<=_order; i++) {
      for (j=0; j <= _order-i; j++) {
          k = _order - i - j;
          int index = TriInterpCoef::ijk_to_index(i,j,k);
          tmp += coef->coef(i,j,k)*_vertices[index]->point();
          x_s += coef->sderiv(i,j,k)*_vertices[index]->point();
          x_r += coef->rderiv(i,j,k)*_vertices[index]->point();
        }
    }
  result->set_point(tmp);

  if (result->has_normal()) {
      for (i=0; i<_order; i++) {
          for (j=0; j <= _order-i; j++) {
              k = _order - i - j;
              int index = TriInterpCoef::ijk_to_index(i,j,k);
              tmp += coef->coef(i,j,k)*_vertices[index]->point();
            }
        }
      result->set_normal(tmp);
    }

  // Find the surface element
  dA = x_r.cross(x_s);
}

void
Triangle::interpolate(double r, double s,
                      const Ref<Vertex>&result, SCVector3&dA,
                      const Ref<Volume> &vol, double isoval)
{
  // set up an initial dummy norm
  SCVector3 norm(0.0);
  result->set_normal(norm);

  // initial guess
  interpolate(r,s,result,dA);

  // now refine that guess
  SCVector3 trialpoint = result->point();
  SCVector3 trialnorm = result->normal();
  SCVector3 newpoint;
  vol->solve(trialpoint,trialnorm,isoval,newpoint);
  // compute the true normal
  vol->set_x(newpoint);
  if (vol->gradient_implemented()) {
      vol->get_gradient(trialnorm);
    }
  trialnorm.normalize();
  result->set_point(newpoint);
  result->set_normal(trialnorm);
}

void
Triangle::flip()
{
  _orientation0 = _orientation0?0:1;
  _orientation1 = _orientation1?0:1;
  _orientation2 = _orientation2?0:1;
}

void
Triangle::set_order(int order, const Ref<Volume>&vol, double isovalue)
{
  if (order > max_order) {
      ExEnv::errn() << scprintf("Triangle::set_order: max_order = %d but order = %d\n",
                       max_order, order);
      abort();
    }

  unsigned int i, j, k;

  if (edge(0)->order() != order
      ||edge(1)->order() != order
      ||edge(2)->order() != order) {
      ExEnv::errn() << "Triangle::set_order: edge order doesn't match" << endl;
      abort();
    }

  _order = order;