ProblemImplicit.cc

#include "ProblemImplicit.h"

namespace AMDiS {
  namespace calcspace {
    void readDofvec(std::istream& in, DOFVector<double>* vec, Mesh* mesh) 
    {
      long size = mesh->getNumberOfVertices();
      std::string buf;
      getline(in, buf);
      while (buf != "vertex values: solution") 
	getline(in, buf);
      for (long i = 0; i < size; i++) {
        in >> buf;
        (*vec)[i] = atof(buf.c_str());
      }
    }


    void readR(std::istream& inStream, double eps, 
	       DOFVector<double>* r, DOFVector<double>* phi1, 
	       DOFVector<double>* phi2, DOFVector<double>* levelSet, 
	       Mesh* mesh) 
    {
      readDofvec(inStream, r, mesh);
      std::vector<double>& vecR = r->getVector();
      std::vector<double>& vecPhi1 = phi1->getVector();
      std::vector<double>& vecPhi2 = phi2->getVector();
      std::vector<double>& vecLevelSet = levelSet->getVector();
      bool checkSize = 
	vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size();
      TEST_EXIT(checkSize)("something went wrong\n");
      for (int i = vecR.size() - 1; i >= 0; --i) {
        vecPhi1[i] = Phi1(vecR[i], eps);
        vecPhi2[i] = Phi2(vecR[i], eps);
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }


    void readPhi1(istream& inStream, double eps, 
		  DOFVector<double>* r, DOFVector<double>* phi1, 
		  DOFVector<double>* phi2, DOFVector<double>* levelSet, 
		  Mesh* mesh) 
    {
      readDofvec(inStream, phi1, mesh);
      std::vector<double>& vecR = r->getVector();
      std::vector<double>& vecPhi1 = phi1->getVector();
      std::vector<double>& vecPhi2 = phi2->getVector();
      std::vector<double>& vecLevelSet = levelSet->getVector();
      
      bool checkSize = 
	vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size();
      TEST_EXIT(checkSize)("something went wrong\n");
      for (int i = vecR.size() - 1; i>= 0; --i) {
        vecR[i] = Phi1ToR(vecPhi1[i], eps);
        //vecPhi2[i] = Phi2(vecR[i], eps);
        vecPhi2[i] = 1 - vecPhi1[i];
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }


    void readPhi2(istream& inStream, double eps, 
		    DOFVector<double>* r, DOFVector<double>* phi1, 
		    DOFVector<double>* phi2, DOFVector<double>* levelSet, 
		    Mesh* mesh) 
    {
      readDofvec(inStream, phi2, mesh);
      std::vector<double>& vecR = r->getVector();
      std::vector<double>& vecPhi1 = phi1->getVector();
      std::vector<double>& vecPhi2 = phi2->getVector();
      std::vector<double>& vecLevelSet = levelSet->getVector();

      bool checkSize = vecR.size() != vecPhi1.size() && vecPhi2.size() != vecR.size();
      TEST_EXIT(checkSize)("something went wrong\n");
      for (int i = vecR.size() - 1; i >= 0; --i) {
        vecR[i] = Phi2ToR(vecPhi2[i], eps);
        //vecPhi1[i] = Phi1(vecR[i], eps);
        vecPhi1[i] = 1 - vecPhi2[i];
        vecLevelSet[i] = LevelSet(vecR[i]);
      }
    }
  }


  bool ProblemImplicitScal::createImplicitMesh() 
  {
    std::string meshFilename("");
    GET_PARAMETER(0, name + "->implicit mesh->mesh file", &meshFilename);
    //if no file name is given, there's no implicit mesh
    if (meshFilename.length() == 0)
      return false;

    std::string dofFilename("");
    GET_PARAMETER(0, name + "->implicit mesh->dof file", &dofFilename); 
    if (dofFilename.length() == 0)
      return false;

    double eps(-1.0);
    GET_PARAMETER(0, name + "->implicit mesh->eps", "%d", &eps);
    if (eps < 0)
      return false;

    int serType(-1);
    GET_PARAMETER(0, name + "->implicit mesh->type", "%d", &serType);
    if (serType < 0)
      return false;

    TEST_EXIT(mesh != NULL)("the mesh was not created\n");

    std::ifstream inStream(meshFilename.c_str());
    mesh->deserialize(inStream);
    inStream.close();
    //create the fespace with the correct admin

    createFeSpace();

    r = new DOFVector< double >(getFeSpace(), "r");
    phi1 = new DOFVector<double>(getFeSpace(), "phi1");
    phi2 = new DOFVector<double>(getFeSpace(), "phi2");
    levelSet = new DOFVector< double >(getFeSpace(), "levelSet");
    inStream.open(dofFilename.c_str());
    switch (serType) {
    case 0:
      calcspace::readR(inStream, eps, r, phi1, phi2, levelSet, mesh);
      break;
    case 1:
      calcspace::readPhi1(inStream, eps, r, phi1, phi2, levelSet, mesh);
      break;
    case 2:
      calcspace::readPhi2(inStream, eps, r, phi1, phi2, levelSet, mesh);
      break;
    default:
      WARNING("unknown implicit mesh type\n");
    }
    inStream.close();
    return true;
  }


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


  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 (int i = 0; i < meshes.size(); ++i)
      createImplicitMesh(i);    
  }


  bool ProblemImplicitVec::createImplicitMesh(int p) 
  {
    std::string path = 
      name + "->implicit mesh[" + boost::lexical_cast< std::string >(p) + "]->";
    std::string meshFilename("");
    GET_PARAMETER(0, path + "mesh file", &meshFilename);
    //if no file name is given, there's no implicit mesh
    if (meshFilename.length() == 0)
      return false;
    std::string dofFilename("");
    GET_PARAMETER(0,path + "dof file", &dofFilename); 
    if (dofFilename.length() == 0)
      return false;

    double eps(-1.0);
    GET_PARAMETER(0, path + "eps", "%d", &eps);
    if (eps < 0.0)
      return false;
    int serType(-1);
    GET_PARAMETER(0, path + "type", "%d", &serType);
    if (serType < 0)
      return false;
    TEST_EXIT(meshes[p] != NULL)("the mesh was not created\n");
    std::ifstream inStream(meshFilename.c_str());
    meshes[p]->deserialize(inStream);
    inStream.close();
    //create the fespace with the correct admin
    createFeSpace(NULL);
    r[p] = new DOFVector< double >(getFeSpace(p), "r");
    phi1[p] = new DOFVector<double>(getFeSpace(p), "phi1");
    phi2[p] = new DOFVector<double>(getFeSpace(p), "phi2");
    levelSet[p] = new DOFVector< double >(getFeSpace(p), "levelSet");
    inStream.open(dofFilename.c_str());
    switch (serType) {
    case 0:
      calcspace::readR(inStream, eps,r[p], phi1[p], phi2[p],
		       levelSet[p], meshes[p]);
      break;
    case 1:
      calcspace::readPhi1(inStream, eps, r[p], phi1[p], phi2[p],
			  levelSet[p], meshes[p]);
      break;
    case 2:
      calcspace::readPhi2(inStream,eps, r[p], phi1[p], phi2[p],
			  levelSet[p], meshes[p]);
      break;
    default:
      WARNING("unknown implicit mesh type\n");
    }
    inStream.close();
    return true;
  }


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