ConformTransformation.cpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane Del Pino

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program 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 General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: ConformTransformation.cpp,v 1.17 2007/06/09 10:37:07 delpinux Exp $

#include <Domain.hpp>

#include <ConformTransformation.hpp>
#include <ScalarFunctionBase.hpp>

real_t
ConformTransformationQ1Hexahedron::integrate(const ScalarFunctionBase& f) const
{
  throw ErrorHandler(__FILE__,__LINE__,
                     "not implemented",
                     ErrorHandler::unexpected);
  return 0;
}

bool
ConformTransformationQ1Hexahedron::
__inside(const real_t& x,
         const real_t& y,
         const real_t& z,
         TinyVector<6, bool>& faces) const
{
  bool inside = true;
  faces = false;
  for (size_t i = 0 ; i<Hexahedron::NumberOfFaces; ++i) {
    TinyVector<3, real_t> faceMassCenter(0,0,0);
    TinyVector<Hexahedron::FaceType::NumberOfVertices, Vertex*> face;

    for (size_t n=0; n<Hexahedron::FaceType::NumberOfVertices; ++n) {
      const Vertex& v = __H(Hexahedron::faces[i][n]);
      faceMassCenter += v;
      face[n] = const_cast<Vertex*>(&v);
    }
    Quadrangle Q(face);
    faceMassCenter *= (1./Hexahedron::FaceType::NumberOfVertices);
    TinyVector<3,real_t> normal = Q.normal();

    TinyVector<3, real_t> X(x,y,z);
    X -= faceMassCenter;
    if (X*normal > 0) { // must look in that direction
      faces[i] = true;
      inside = false;
    }
 }
  return inside;
}

bool ConformTransformationQ1Hexahedron::
invertT(const real_t& x,
        const real_t& y,
        const real_t& z,
        TinyVector<3,real_t>& Xhat) const
{
  // initialization
  for(size_t i=0; i<3; ++i)
    Xhat[i] = 0.5;

  TinyVector<3,real_t> X;
  TinyVector<3,real_t> f;
  TinyVector<3,real_t> delta;
  const size_t maxiter = 100;
  size_t niter = 0;

  do {
    niter++;

    value(Xhat,f);
    f[0] -= x;
    f[1] -= y;
    f[2] -= z;

    // Evaluation of the Jacobian 
    TinyMatrix<3,3,real_t> J;
    dx(Xhat,X);
    for (size_t i=0; i<3; ++i)
      J(i,0) = X[i];

    dy(Xhat,X);
    for (size_t i=0; i<3; ++i)
      J(i,1) = X[i];

    dz(Xhat,X);
    for (size_t i=0; i<3; ++i)
      J(i,2) = X[i];

    delta = f/J;

    // Uses relaxation to help convergence (the functional is not convexe)
    Xhat -= delta;
    
    if (niter>maxiter) { // or(Norm(f)>10)) {
      return false;
    }

  } while(Norm(f)>1E-3);

  for (size_t i = 0; i<Xhat.size(); ++i) {
    if ((Xhat[i]<1E-3) or (Xhat[i]>1.001)) {
      return false;
    }
  }

  return true;
}



real_t
ConformTransformationP1Tetrahedron::integrate(const ScalarFunctionBase& f) const
{
  throw ErrorHandler(__FILE__,__LINE__,
                     "not implemented",
                     ErrorHandler::unexpected);
  return 0;
}

real_t
ConformTransformationQ1CartesianHexahedron::
integrateCharacteristic(const Domain& d) const
{
  return 1;
  // Here we use order 4 Lobatto quadrature

  TinyVector<4, real_t>  x;
  x[0] = 0.;
  x[1] = .27639320225002103036; //(1-sqrt(5)/5)/2.;
  x[2] = .72360679774997896964; //(1+sqrt(5)/5)/2.;
  x[3] = 1.;

  TinyVector<4, real_t>  w;
  w[0] = 1./12.;
  w[1] = 5./12.;
  w[2] = 5./12.;
  w[3] = 1./12.;

  real_t sum = 0;
  TinyVector<3, real_t> X_hat;
  TinyVector<3, real_t> X;
  for (unsigned i=0; i<4; ++i) {
    X_hat[0] = x[i];
    for (unsigned j=0; j<4; ++j) {
      X_hat[1] = x[j];
      for (unsigned k=0; k<4; ++k) {
        X_hat[2] = x[k];
        this->value(X_hat, X);
        sum += w[i]*w[j]*w[k] * (d.inside(X) ? 1 : 0);
      }
    }
  }
  return sum;
}


real_t
ConformTransformationP1Triangle::integrate(const ScalarFunctionBase& f) const
{
  throw ErrorHandler(__FILE__,__LINE__,
                     "not implemented",
                     ErrorHandler::unexpected);

  return 0.;
}

real_t
ConformTransformationQ1Quadrangle::integrate(const ScalarFunctionBase& f) const
{
  throw ErrorHandler(__FILE__,__LINE__,
                     "not implemented",
                     ErrorHandler::unexpected);
  return 0.;
}

---