ProblemImplicit.cc

#include "ProblemImplicit.h"
#include "MathFunctions.h"

using namespace std;
namespace AMDiS {
  void ProblemImplicitBase::readDofVec(std::istream& in, DOFVector<double>* vec, 
                                       Mesh* mesh) 
  {
    long size = mesh->getNumberOfVertices();

    TEST_EXIT(vec->getSize()>=size)("dofvector smaller than vertex data vector\n");
    std::string buf;
    getline(in, buf);
    while (!in.eof() && buf != "vertex values: solution") 
      getline(in, buf);

    TEST_EXIT(!in.eof())("no vertex data\n");

    for (long i = 0; i < size ; ++i) {
      in >> buf;
      (*vec)[i] = atof(buf.c_str());
    }
  }

  void ProblemImplicitBase::readR(std::istream& inStream, double eps, 
	                          Mesh* mesh, int implMesh , int comp) 
  {
    FUNCNAME("ProblemImplicitBase::readR()");

    DOFVector< double >* r = getSignedDistance(implMesh, comp);
    DOFVector< double >* phi1 = getPhi1(implMesh, comp);
    DOFVector< double >* phi2 = getPhi2(implMesh, comp);
    DOFVector< double >* levelSet = getLevelset(implMesh, comp);

    TEST_EXIT( r != NULL )("no signed distance vector\n");
    TEST_EXIT( phi1 != NULL )("no phasefield1 vector\n");
    TEST_EXIT( phi2 != NULL )("no phasefield2 vector\n");
    TEST_EXIT( levelSet != NULL )("no levelSet vector\n");

    bool checkSize = r->getSize() == phi1->getSize() && 
	             r->getSize() == phi2->getSize();

    TEST_EXIT(checkSize)("something went wrong\n");
    

    readDofVec(inStream, r, mesh);

    for (int i = r->getSize() - 1; i >= 0; --i) {
      (*phi1)[i] = Phi1((*r)[i], eps);
      (*phi2)[i] = Phi2((*r)[i], eps);
      (*levelSet)[i] = LevelSet((*r)[i]);
    }
  }


  void ProblemImplicitBase::readPhi1(istream& inStream, double eps, 
                                     Mesh* mesh, int implMesh, int comp )
  {
    DOFVector< double >* r = getSignedDistance(implMesh, comp);
    DOFVector< double >* phi1 = getPhi1(implMesh, comp);
    DOFVector< double >* phi2 = getPhi2(implMesh, comp);
    DOFVector< double >* levelSet = getLevelset(implMesh, comp);

    TEST_EXIT( r != NULL )("no signed distance vector\n");
    TEST_EXIT( phi1 != NULL )("no phasefield1 vector\n");
    TEST_EXIT( phi2 != NULL )("no phasefield2 vector\n");
    TEST_EXIT( levelSet != NULL )("no levelSet vector\n");

    bool checkSize = r->getSize() == phi1->getSize() && 
	             r->getSize() == phi2->getSize();

    TEST_EXIT(checkSize)("something went wrong\n");

    readDofVec(inStream, phi1, mesh);

    for (int i = r->getSize() - 1; i>= 0; --i) {
      (*r)[i] = Phi1ToR((*phi1)[i], eps);
      (*phi2)[i] = 1 - (*phi1)[i];
      (*levelSet)[i] = LevelSet((*r)[i]);
    }
  }


  void ProblemImplicitBase::readPhi2(istream& inStream, double eps, 
	                             Mesh* mesh, int implMesh, int comp )
  { 
    DOFVector< double >* r = getSignedDistance(implMesh, comp);
    DOFVector< double >* phi1 = getPhi1(implMesh, comp);
    DOFVector< double >* phi2 = getPhi2(implMesh, comp);
    DOFVector< double >* levelSet = getLevelset(implMesh, comp);

    TEST_EXIT( r != NULL )("no signed distance vector\n");
    TEST_EXIT( phi1 != NULL )("no phasefield1 vector\n");
    TEST_EXIT( phi2 != NULL )("no phasefield2 vector\n");
    TEST_EXIT( levelSet != NULL )("no levelSet vector\n");

    bool checkSize = r->getSize() == phi1->getSize() &&
	             r->getSize() == phi2->getSize();
    TEST_EXIT(checkSize)("something went wrong\n");

    readDofVec(inStream, phi2, mesh);

    for (int i = r->getSize() - 1; i >= 0; --i) {
      (*r)[i] = Phi2ToR((*r)[i], eps);
      (*phi1)[i] = 1 - (*phi2)[i];
      (*levelSet)[i] = LevelSet((*r)[i]);
    }
  }
  
  std::string ProblemImplicitBase::getDofFilename(std::string path, int implMesh)
  {
    std::string suffix = "[" + 
	                 boost::lexical_cast< std::string >(implMesh) +
			 "]";

    std::string dofFilename("");
    GET_PARAMETER(0, path + "dof file" + suffix, &dofFilename);
    if (implMesh == 0 && dofFilename.length() == 0)
      GET_PARAMETER(0, path + "dof file", &dofFilename);
    return dofFilename;    
  }

  double ProblemImplicitBase::getEpsilon(std::string path, int implMesh)
  {
    std::string suffix = "[" + 
	                 boost::lexical_cast< std::string >(implMesh) +
			 "]";

    double eps(-1);
    GET_PARAMETER(0, path + "eps" + suffix, "%g", &eps);
    if (implMesh == 0 && eps < 0)
      GET_PARAMETER(0, path + "eps", "%g", &eps);
    return eps;
  }

  int ProblemImplicitBase::getType(std::string path, int implMesh)
  {
    std::string suffix = "[" + 
	                 boost::lexical_cast< std::string >(implMesh) +
			 "]";
    int serType(-1);
    GET_PARAMETER(0, path + "type" + suffix, "%d", &serType);
    if (implMesh == 0 && serType < 0)
      GET_PARAMETER(0, path + "type", "%d", &serType);
    return serType;
  }

  DOFVector< double >* ProblemImplicitScal::getSignedDistance(int im , int )
  { 
    if (readImplMesh)
      return r[im]; 
    return NULL;
  }

  DOFVector< double >* ProblemImplicitScal::getPhi1(int im , int )
  { 
    if (readImplMesh)
      return phi1[im];
    return NULL;
  }

  DOFVector< double >* ProblemImplicitScal::getPhi2(int im , int )
  { 
    if (readImplMesh)
      return phi2[im]; 
    return NULL;
  }

  DOFVector< double >* ProblemImplicitScal::getLevelset(int im , int ) 
  { 
    if (readImplMesh)
      return levelSet[im]; 
    return NULL;
  }

  void ProblemImplicitScal::createMesh()
  {
    FUNCNAME("ProblemImplicitScal::createMesh()");
    ProblemScal::createMesh();
    std::string path = name + "->implicit mesh";
    std::string meshFilename("");

    GET_PARAMETER(0, path + "->macro file name", &meshFilename);
    std::string serFilename("");
    GET_PARAMETER(0, path + "->serialization file name", &serFilename);
    if ( meshFilename.length() > 0) {
      mesh->setName(path);
      readImplMesh = true;
    }
    else 
      if ( serFilename.length() > 0 ) {
        //INFO(5,1)("loading implicit mesh %s \n", serFilename.c_str());
	std::string oldName = mesh->getName();
        std::ifstream inStream(serFilename.c_str());
        mesh->deserialize(inStream);
        inStream.close();
	mesh->setName(oldName);
        readImplMesh = true;
      }

  }

  void ProblemImplicitScal::initialize(Flag initFlag, ProblemScal* adaptProblem, 
		                       Flag adoptFlag)
  {
     FUNCNAME("ProblemImplicitScal::initialize()");
     ProblemScal::initialize(initFlag);
     createImplicitMesh();
  }

  bool ProblemImplicitScal::createImplicitMesh() 
  {
    FUNCNAME("ProblemImplicitScal::createImplicitMesh()");
    INFO(5,1)("setting implicit mesh\n");
    if (!isImplicitMesh())
      return false;
    std::string path = name + "->implicit mesh->";
    int nImplMesh(1);
    GET_PARAMETER(0, name + "nr meshes", "%d", &nImplMesh);

    r.resize(nImplMesh);
    phi1.resize(nImplMesh);
    phi2.resize(nImplMesh);
    levelSet.resize(nImplMesh);
    for ( int i = 0; i < nImplMesh ; ++i) {
      r[i] = new DOFVector< double >(getFeSpace(), "r");
      phi1[i] = new DOFVector< double >(getFeSpace(), "phi1");
      phi2[i] = new DOFVector< double >(getFeSpace(), "phi2");
      levelSet[i] = new DOFVector< double >(getFeSpace(), "levelSet");
      createImplicitMesh(path, i);
    }
    return true;
  }

  bool ProblemImplicitScal::createImplicitMesh(string path, int p)
  {
    
    std::string dofFilename = getDofFilename(path, p);
    if (dofFilename.length() == 0)
      return false;

    double eps = getEpsilon(path, p);
    if (eps < 0.0)
      return false;
    int serType = getType(path, p);
    if (serType < 0)
      return false;

    TEST_EXIT(mesh != NULL)("the mesh was not created\n");
    
    std::ifstream inStream(dofFilename.c_str());
    switch (serType) {
    case 0:
      readR(inStream, eps, mesh, p);
      break;
    case 1:
      readPhi1(inStream, eps, mesh, p);
      break;
    case 2:
      readPhi2(inStream, eps, mesh, p);
      break;
    default:
      WARNING("unknown implicit mesh type\n");
    }
    inStream.close();
    return true;
  }


  DOFVector< double >* ProblemImplicitVec::getSignedDistance(unsigned int im , unsigned int m) 
  { 
    if ( m >= r.size() || im >= r[m].size())
	return NULL;
      return (r[m])[im]; 
    }

  DOFVector< double >* ProblemImplicitVec::getPhi1(unsigned int im, unsigned int m)
  {
    if (m >= phi1.size() || im >= phi1[m].size())
      return NULL;

    return (phi1[m])[im];
  }

  DOFVector< double >* ProblemImplicitVec::getPhi2(unsigned int im, unsigned int m)
  {
    if (m >= phi2.size() || im >= phi2[m].size())
      return NULL;

    return (phi2[m])[im];
  }

  DOFVector< double >* ProblemImplicitVec::getLevelset(unsigned int im, unsigned int m)
  {
    if (m >= levelSet.size() || im >= levelSet[m].size())
      return NULL;

    return (levelSet[m])[im];
  }

  bool ProblemImplicitVec::createImplicitMesh()
  {
    //check each mesh if it's an implicit one
    r.resize(meshes.size());
    phi1.resize(meshes.size());
    phi2.resize(meshes.size());
    levelSet.resize(meshes.size());
    for (unsigned int i = 0; i < meshes.size(); ++i)
      createImplicitMesh(i);
    return true;
  }


  bool ProblemImplicitVec::createImplicitMesh(int p) 
  {
    FUNCNAME("ProblemImplicitVec::createImplicitMesh()");
    //INFO(5,1)("initing mesh %d\n", p);
    if (!isImplicitMesh(p))
      return false;
    std::string path = name + "->implicit mesh[" 
	             + boost::lexical_cast< std::string >(p) + "]->";
    int nImplMeshes(0);
    GET_PARAMETER(0, path + "nr meshes", "%d", &nImplMeshes);
    //INFO(5,1)("has %d meshe(s)\n", nImplMeshes);
    r[p].resize(nImplMeshes);
    phi1[p].resize(nImplMeshes);
    phi2[p].resize(nImplMeshes);
    levelSet[p].resize(nImplMeshes);

    for ( int i = 0; i < nImplMeshes ; ++i ) {
      (r[p])[i] = new DOFVector< double >(getFeSpace(p), "r");
      (phi1[p])[i] = new DOFVector< double >(getFeSpace(p), "phi1");
      (phi2[p])[i] = new DOFVector< double >(getFeSpace(p), "phi2");
      (levelSet[p])[i] = new DOFVector< double >(getFeSpace(p), "levelSet");
      createImplicitMesh(path, i, p);
    }
    return true;
  }
  
  bool ProblemImplicitVec::createImplicitMesh(std::string path, int implMesh, 
		                              int comp)
  {
    FUNCNAME("ProblemImplicitVec::createImplicitMesh()");
    std::string dofFilename = getDofFilename(path, implMesh);
    if (dofFilename.length() == 0)
      return false;

    double eps = getEpsilon(path, implMesh);
    if (eps < 0.0)
      return false;
    int serType = getType(path, implMesh);
    if (serType < 0)
      return false;

    TEST_EXIT(meshes[comp] != NULL)("the mesh was not created\n");
    
    std::ifstream inStream(dofFilename.c_str());
    switch (serType) {
    case 0:
      readR(inStream, eps, meshes[comp], implMesh, comp);
      break;
    case 1:
      readPhi1(inStream, eps, meshes[comp], implMesh, comp);
      break;
    case 2:
      readPhi2(inStream, eps, meshes[comp], implMesh, comp);
      break;
    default:
      WARNING("unknown implicit mesh type\n");
    }
    inStream.close();
    return true;
  }

  void ProblemImplicitVec::createMesh()
  {
    FUNCNAME("ProblemImplicitVec::createMesh()");
    ProblemVec::createMesh();
    implMesh.resize(meshes.size());
    for (unsigned int i = 0 ; i < meshes.size() ; ++i )
    {
      std::string path = name + "->implicit mesh[" +
                         boost::lexical_cast< std::string >(i) + "]";
      std::string meshFilename("");
      GET_PARAMETER(0, path + "->macro file name", &meshFilename);
      std::string serFilename("");
      GET_PARAMETER(0, path + "->serialization file name", &serFilename);
      implMesh[i] = true;
      if (meshFilename.length() > 0)
        meshes[i]->setName(path);
      else if (serFilename.length() > 0) 
      {
        //INFO(5,1)("reading serialization file name %s \n", serFilename.c_str());
	std::ifstream inStream(serFilename.c_str());
	meshes[i]->deserialize(inStream);
	inStream.close();
      } else
        implMesh[i] = false;
    }

  }

  void ProblemImplicitVec::initialize(Flag initFlag, ProblemScal* adaptProblem, 
		                      Flag adoptFlag)
  {
     FUNCNAME("ProblemImplicitVec::initialize()");
     ProblemVec::initialize(initFlag);    
     createImplicitMesh();
  }
}