triangle.cc

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

  delete[] _vertices;

  _vertices = new Ref<Vertex>[TriInterpCoef::order_to_nvertex(_order)];

  // fill in the corner vertices
  _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);

  // fill in the interior edge vertices
  for (i = 1; i < _order; i++) {
      j = _order - i;
      _vertices[TriInterpCoef::ijk_to_index(0, i, j)]
          = _edges[1]->interior_vertex(_orientation1?i:j);
      _vertices[TriInterpCoef::ijk_to_index(j, 0, i)]
          = _edges[0]->interior_vertex(_orientation0?i:j);
      _vertices[TriInterpCoef::ijk_to_index(i, j, 0)]
          = _edges[2]->interior_vertex(_orientation2?i:j);
    }

  const SCVector3& p0 = vertex(0)->point();
  const SCVector3& p1 = vertex(1)->point();
  const SCVector3& p2 = vertex(2)->point();
  const SCVector3& norm0 = vertex(0)->normal();
  const SCVector3& norm1 = vertex(1)->normal();
  const SCVector3& norm2 = vertex(2)->normal();

  for (i=0; i<=_order; i++) {
      double I = (1.0*i)/_order;
      for (j=0; j<=_order-i; j++) {
          SCVector3 trialpoint;
          SCVector3 trialnorm;
          SCVector3 newpoint;
          double J = (1.0*j)/_order;
          k = _order - i - j;
          if (!i || !j || !k) continue; // interior point check
          double K = (1.0*k)/_order;
          int index = TriInterpCoef::ijk_to_index(i,j,k);
          // this get approximate vertices and normals
          trialpoint = I*p0 + J*p1 + K*p2;
          trialnorm = I*norm0 + J*norm1 + K*norm2;
          // now refine that guess
          vol->solve(trialpoint,trialnorm,isovalue,newpoint);
          // compute the true normal
          vol->set_x(newpoint);
          if (vol->gradient_implemented()) {
              vol->get_gradient(trialnorm);
            }
          trialnorm.normalize();
          _vertices[index] = new Vertex(newpoint,trialnorm);
        }
    }
}

/////////////////////////////////////////////////////////////////////////
// TriangleIntegrator

static ClassDesc TriangleIntegrator_cd(
  typeid(TriangleIntegrator),"TriangleIntegrator",1,"public DescribedClass",
  0, create<TriangleIntegrator>, 0);

TriangleIntegrator::TriangleIntegrator(int order):
  _n(order)
{
  _r = new double [_n];
  _s = new double [_n];
  _w = new double [_n];
  coef_ = 0;
}

TriangleIntegrator::TriangleIntegrator(const Ref<KeyVal>& keyval)
{
  _n = keyval->intvalue("n");
  if (keyval->error() != KeyVal::OK) {
      _n = 7;
    }
  _r = new double [_n];
  _s = new double [_n];
  _w = new double [_n];
  coef_ = 0;
}

TriangleIntegrator::~TriangleIntegrator()
{
  delete[] _r;
  delete[] _s;
  delete[] _w;
  clear_coef();
}

void
TriangleIntegrator::set_n(int n)
{
  delete[] _r;
  delete[] _s;
  delete[] _w;
  _n = n;
  _r = new double [_n];
  _s = new double [_n];
  _w = new double [_n];
  clear_coef();
}

void
TriangleIntegrator::set_w(int i,double w)
{
  _w[i] = w;
}

void
TriangleIntegrator::set_r(int i,double r)
{
  _r[i] = r;
}

void
TriangleIntegrator::set_s(int i,double s)
{
  _s[i] = s;
}

void
TriangleIntegrator::init_coef()
{
  int i, j;

  clear_coef();

  coef_ = new Ref<TriInterpCoef>*[Triangle::max_order];
  for (i=0; i<Triangle::max_order; i++) {
      coef_[i] = new Ref<TriInterpCoef>[_n];
      for (j=0; j<_n; j++) {
          TriInterpCoefKey key(i+1,_r[j],_s[j]);
          coef_[i][j] = new TriInterpCoef(key);
        }
    }
}

void
TriangleIntegrator::clear_coef()
{
  int i, j;

  if (coef_) {
      for (i=0; i<Triangle::max_order; i++) {
          for (j=0; j<_n; j++) {
              coef_[i][j] = 0;
            }
          delete[] coef_[i];
        }
    }
  delete[] coef_;
  coef_ = 0;
}

/////////////////////////////////////////////////////////////////////////
// GaussTriangleIntegrator

static ClassDesc GaussTriangleIntegrator_cd(
  typeid(GaussTriangleIntegrator),"GaussTriangleIntegrator",1,"public TriangleIntegrator",
  0, create<GaussTriangleIntegrator>, 0);

GaussTriangleIntegrator::GaussTriangleIntegrator(const Ref<KeyVal>& keyval):
  TriangleIntegrator(keyval)
{
  ExEnv::out0() << "Created a GaussTriangleIntegrator with n = " << n() << endl;
  init_rw(n());
  init_coef();
}

GaussTriangleIntegrator::GaussTriangleIntegrator(int order):
  TriangleIntegrator(order)
{
  init_rw(n());
  init_coef();
}

void
GaussTriangleIntegrator::set_n(int n)
{
  TriangleIntegrator::set_n(n);
  init_rw(n);
  init_coef();
}

void
GaussTriangleIntegrator::init_rw(int order)
{
  if (order == 1) {
      set_r(0, 1.0/3.0);
      set_s(0, 1.0/3.0);
      set_w(0, 1.0);
    }
  else if (order == 3) {
      set_r(0, 1.0/6.0);
      set_r(1, 2.0/3.0);
      set_r(2, 1.0/6.0);
      set_s(0, 1.0/6.0);
      set_s(1, 1.0/6.0);
      set_s(2, 2.0/3.0);
      set_w(0, 1.0/3.0);
      set_w(1, 1.0/3.0);
      set_w(2, 1.0/3.0);
    }
  else if (order == 4) {
      set_r(0, 1.0/3.0);
      set_r(1, 1.0/5.0);
      set_r(2, 3.0/5.0);
      set_r(3, 1.0/5.0);
      set_s(0, 1.0/3.0);
      set_s(1, 1.0/5.0);
      set_s(2, 1.0/5.0);
      set_s(3, 3.0/5.0);
      set_w(0, -27.0/48.0);
      set_w(1, 25.0/48.0);
      set_w(2, 25.0/48.0);
      set_w(3, 25.0/48.0);
    }
  else if (order == 6) {
      set_r(0, 0.091576213509771);
      set_r(1, 0.816847572980459);
      set_r(2, r(0));
      set_r(3, 0.445948490915965);
      set_r(4, 0.108103018168070);
      set_r(5, r(3));
      set_s(0, r(0));
      set_s(1, r(0));
      set_s(2, r(1));
      set_s(3, r(3));
      set_s(4, r(3));
      set_s(5, r(4));
      set_w(0, 0.109951743655322);
      set_w(1, w(0));
      set_w(2, w(0));
      set_w(3, 0.223381589678011);
      set_w(4, w(3));
      set_w(5, w(3));
    }
  else if (order == 7) {
      set_r(0, 1.0/3.0);
      set_r(1, 0.101286507323456);
      set_r(2, 0.797426985353087);
      set_r(3, r(1));
      set_r(4, 0.470142064105115);
      set_r(5, 0.059715871789770);
      set_r(6, r(4));
      set_s(0, r(0));
      set_s(1, r(1));
      set_s(2, r(1));
      set_s(3, r(2));
      set_s(4, r(4));
      set_s(5, r(4));
      set_s(6, r(5));
      set_w(0, 0.225);
      set_w(1, 0.125939180544827);
      set_w(2, w(1));
      set_w(3, w(1));
      set_w(4, 0.132394152788506);
      set_w(5, w(4));
      set_w(6, w(4));
    }
  else {
      ExEnv::errn() << "GaussTriangleIntegrator: invalid order " << order << endl;
      abort();
    }

  // scale the weights by the area of the unit triangle
  for (int i=0; i<order; i++) {
      set_w(i, w(i)*0.5);
    }
}

GaussTriangleIntegrator::~GaussTriangleIntegrator()
{
}

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

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