amroc/amr/BrepOutput.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2009 Ralf Deiterding

#ifndef AMROC_BREP_OUTPUT_H
#define AMROC_BREP_OUTPUT_H

#include "AMRInterpolation.h"
#include "geom/kernel.h"

template <class VectorType, int dim> class BrepOutOps;

template <class VectorType>
class BrepOutOps<VectorType,2> {
public:
  BrepOutOps() {}

  void ElementNormals(const int num_pos, const ads::FixedArray<2,DataType>* pos,
                      const int num_elem, const ads::FixedArray<2,int>* conn,
                      ads::FixedArray<2,DataType>* normals) {
    for (register int i=0; i<num_elem; i++) {
      ads::FixedArray<2,DataType> tangent(pos[conn[i][1]]);
      tangent -= pos[conn[i][0]];
      geom::normalize(&tangent);
      normals[i][0] = -tangent[1];
      normals[i][1] =  tangent[0];
    }
  }

  void ElementAreas(const int num_pos, const ads::FixedArray<2,DataType>* pos,
                    const int num_elem, const ads::FixedArray<2,int>* conn,
                    DataType* areas) {
    for (register int i=0; i<num_elem; i++)
      areas[i] = geom::computeDistance(pos[conn[i][0]],pos[conn[i][1]]); 
  }

  void CrossProducts(const int num_vec, const ads::FixedArray<2,DataType>* v1, 
                     const ads::FixedArray<2,DataType>* v2, 
                     ads::FixedArray<2,DataType>* r) {
    for (register int i=0; i<num_vec; i++) {
      r[i][0] = v1[i][0]*v2[i][1]-v1[i][1]*v2[i][0]; r[i][1] = 0.;
    }
  }

  void DotProducts(const int num_vec, const ads::FixedArray<2,DataType>* v1, 
                   const ads::FixedArray<2,DataType>* v2, 
                   DataType* r) {
    for (register int i=0; i<num_vec; i++) 
      r[i] = geom::computeDotProduct(v1[i],v2[i]);
  }
};


template <class VectorType>
class BrepOutOps<VectorType,3> {
public:
  BrepOutOps() {}

  void ElementNormals(const int num_pos, const ads::FixedArray<3,DataType>* pos,
                      const int num_elem, const ads::FixedArray<3,int>* conn,
                      ads::FixedArray<3,DataType>* normals) {
    for (register int i=0; i<num_elem; i++) {
      const ads::FixedArray<3,int>& face = conn[i];
      geom::computeCrossProduct(pos[face[2]] - pos[face[1]],
                                pos[face[1]] - pos[face[0]], &(normals[i]));
      geom::normalize(&(normals[i]));
    }
  }

  void ElementAreas(const int num_pos, const ads::FixedArray<3,DataType>* pos,
                    const int num_elem, const ads::FixedArray<3,int>* conn,
                    DataType* areas) {
    for (register int i=0; i<num_elem; i++)
      areas[i] = geom::computeArea(pos[conn[i][0]],pos[conn[i][1]],pos[conn[i][2]]); 
  }

  void CrossProducts(const int num_vec, const ads::FixedArray<3,DataType>* v1, 
                     const ads::FixedArray<3,DataType>* v2, 
                     ads::FixedArray<3,DataType>* r) {
    for (register int i=0; i<num_vec; i++) 
      r[i] = geom::computeCrossProduct(v1[i],v2[i]);
  }

  void DotProducts(const int num_vec, const ads::FixedArray<3,DataType>* v1, 
                   const ads::FixedArray<3,DataType>* v2, 
                   DataType* r) {
    for (register int i=0; i<num_vec; i++) 
      r[i] = geom::computeDotProduct(v1[i],v2[i]);
  }
};


template <class VectorType, int dim>
class BrepOutput : public AMRBase<VectorType,dim>, public BrepOutOps<VectorType,dim> {
  typedef AMRBase<VectorType,dim> base;
  typedef typename VectorType::InternalDataType DataType;

public:
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef GridFunction<DataType,dim> grid_fct_type;
  typedef AMRInterpolation<VectorType,dim> interpolation_type;
  typedef typename interpolation_type::point_type point_type;
  typedef ads::FixedArray<dim,DataType> geom_point_type;
  typedef ads::FixedArray<dim,int> multi_index_type;
  typedef BrepOutOps<VectorType,dim> ops;

  BrepOutput(interpolation_type* interpol) : base(), _Interpolation(interpol), 
                                             _Ncnt(0), _OutputType(0), StepFormat(0), 
                                             StepDirectories(0) {
    BaseFileName = "-";
    CompName = (std::string*) 0;
    CompDir = "-";
  } 

  virtual ~BrepOutput() {
    if (CompName) delete [] CompName;
    if (_Interpolation) delete _Interpolation;
  }

  virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    LocCtrl = Ctrl.getSubDevice(prefix+"BrepFile");    
    RegisterAt(LocCtrl,"OutputName",BaseFileName);
    RegisterAt(LocCtrl,"Type",_OutputType);
    if (CompName) delete [] CompName;
    CompName = new std::string[MAXCOMPONENTS];
    for (int d=0; d<MAXCOMPONENTS; d++) {
      char Name[16];
      std::sprintf(Name,"-");
      CompName[d] = Name;
      std::sprintf(Name,"CompName(%d)",d+1);
      RegisterAt(LocCtrl,Name,CompName[d]);
    }    
    RegisterAt(LocCtrl,"OutputDirectory",CompDir);
    RegisterAt(LocCtrl,"StepFormat",StepFormat);
    RegisterAt(LocCtrl,"StepDirectories",StepDirectories);    
  }
  virtual void update() {
    int d=MAXCOMPONENTS;
    for (; d>0; d--) 
      if (CompName[d-1].c_str()[0] != '-') 
        break;
    if (d>0) SetNcnt(d);
  }

  virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
    base::SetupData(gh, ghosts);
    Interpolation().SetupData(gh, ghosts);
  }

  virtual void WriteOut(vec_grid_fct_type& u, 
                        const int num_vertices, const DataType* vertices,
                        const int num_connections, const int* connections) {

    VectorType* data = new VectorType[num_vertices];
    int Time = CurrentTime(base::GH(),0);
    double t = GetPhysicalTime(u,Time,0);
    Interpolation().PointsValues(u,t,num_vertices,(const point_type*)vertices,data);
    int NValues = num_vertices*base::NEquations();
    Interpolation().ArrayCombine(VizServer,NValues,(DataType *)data);
    if (MY_PROC==VizServer) 
      WriteMesh(num_vertices,vertices,num_connections,connections,Time,t,(DataType *)data);
    delete [] data;
  } 

  virtual void WriteOut(grid_fct_type& IOfunc, 
                        const int num_vertices, const DataType* vertices,
                        const int num_connections, const int* connections) {

    DataType* data = new DataType[num_vertices];
    int Time = CurrentTime(base::GH(),0);
    double t = GetPhysicalTime(IOfunc,Time,0);
    Interpolation().PointsValues(IOfunc,t,num_vertices,(const point_type*)vertices,data);
    Interpolation().ArrayCombine(VizServer,num_vertices,data);
    if (MY_PROC==VizServer) 
      WriteMesh(num_vertices,vertices,num_connections,connections,Time,t,data);
    delete [] data;
  } 

  void WriteMesh(const int num_vertices, const DataType* vertices,
                 const int num_connections, const int* connections, 
                 const int Time, const double t, const DataType* data) {

    char ioname[DAGHBktGFNameWidth+256];
    OutputFileName(ioname,BaseFileName.c_str(),Time);
    std::cout << " *** Writing " << ioname << " at t = " << t << std::endl;

    std::ofstream outfile;
    std::ostream* out;
    outfile.open(ioname, std::ios::out);
    out = new std::ostream(outfile.rdbuf());

    if (OutputType()==0) {
      *out << "title = \"mesh t=" << t << "\"" << std::endl;
      *out << "variables = x y z";
      for (int cnt=0;cnt<Ncnt(); cnt++) 
        if (OutputName(cnt).c_str()[0] != '-') 
          *out << " " << OutputName(cnt);
      *out << std::endl;
      *out << "zone t=\"zone 1\", I=" << num_vertices << ", J="
           << num_connections << ", solutiontime=" << t << ", F=FEPOINT, ET=TRIANGLE" << std::endl;
      for (register int j=0;j<num_vertices; j++) {
        if (dim==3) 
          *out << vertices[3*j] << " " << vertices[3*j+1] << " " << vertices[3*j+2]; 
        else if (dim==2)
          *out << vertices[2*j] << " " << vertices[2*j+1] << " 0.";     
        for (int cnt=0;cnt<Ncnt(); cnt++) 
          if (OutputName(cnt).c_str()[0] != '-') 
            *out << " " << data[cnt*num_vertices+j];
        *out << std::endl;
      }
      for (register int j=0;j<num_connections; j++)
        if (dim==3) 
          *out << connections[3*j]+1 << " " << connections[3*j+1]+1 << " " << connections[3*j+2]+1 << std::endl; 
        else if (dim==2)
          *out << connections[2*j]+1 << " " << connections[2*j+1]+1 << " " << connections[2*j+1]+1 << std::endl; 
    }

    outfile.close();
    delete out; 
  }   

  void OutputFileName(char *ioname, const char* name, const int& Time) {
    char time[20], help[DAGHBktGFNameWidth+256];
    if (StepFormat<=0)
      std::sprintf(time,"%d",Time); 
    else if (StepFormat==1) {
      if (Time<=999999) 
        std::sprintf(time,"%6.6d",Time);
      else if (Time<=99999999) 
        std::sprintf(time,"%8.8d",Time);
      else
        std::sprintf(time,"%d",Time);
    }
    if (CompDir.c_str()[0] != '-' && CompDir.length()>0) 
      std::sprintf(ioname,"%s/",CompDir.c_str());
    else
      ioname[0]='\0';
    if (StepDirectories>0) {
      std::sprintf(help,"%s/",time); 
      std::strcat(ioname,help);
    }
    if (OutputType()==0) 
      std::sprintf(help,"%s_%s.tec",name,time); 
    std::strcat(ioname,help);
  }

  DataType Integrate(const DataType* data,
                     const int num_pos, const DataType* vertices,
                     const int num_elem, const int* connections) {
    DataType sum = 0.;
    geom_point_type* area = new DataType[num_elem];
      
    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ops::ElementAreas(num_pos,pos,num_elem,conn,area);
    
    for (register int i=0; i<num_elem; i++)
      if (data[i] != Interpolation().ErrorValue())
        sum += data[i]*area[i];
    
    delete [] area;

    return sum;
  } 

  DataType Integrate(grid_fct_type& IOfunc, 
                     const int num_pos, const DataType* vertices,
                     const int num_elem, const int* connections) {
    geom_point_type* centroids = new geom_point_type[num_elem];
    DataType* data = new DataType[num_elem];

    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ElementCentroids(num_pos,pos,num_elem,conn,centroids);

    int Time = CurrentTime(base::GH(),0);
    double t = GetPhysicalTime(IOfunc,Time,0);
    Interpolation().PointsValues(IOfunc,t,num_elem,reinterpret_cast<const point_type*>(centroids),data);
    Interpolation().ArrayCombine(VizServer,num_elem,data);

    DataType sum = 0.;
    if (MY_PROC==VizServer) {
      geom_point_type* area = new DataType[num_elem];
      
      ops::ElementAreas(num_pos,pos,num_elem,conn,area);
      for (register int i=0; i<num_elem; i++)
        if (data[i] != Interpolation().ErrorValue())
          sum += data[i]*area[i];

     delete [] area;
    }

    delete [] centroids;
    delete [] data;

    return sum;
  } 

  geom_point_type IntegrateNormal(const DataType* data, 
                                  const int num_pos, const DataType* vertices,
                                  const int num_elem, const int* connections) {
    geom_point_type sum(0.);
    DataType* area = new DataType[num_elem];
    geom_point_type* normal = new geom_point_type[num_elem];
    
    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ops::ElementAreas(num_pos,pos,num_elem,conn,area);
    ops::ElementNormals(num_pos,pos,num_elem,conn,normal);
    
    for (register int i=0; i<num_elem; i++) 
      if (data[i] != Interpolation().ErrorValue())
        sum += data[i]*area[i]*normal[i];
    
    delete [] area;
    delete [] normal;

    return sum;
  } 

  geom_point_type IntegrateNormal(grid_fct_type& IOfunc, 
                                  const int num_pos, const DataType* vertices,
                                  const int num_elem, const int* connections) {
    geom_point_type* centroids = new geom_point_type[num_elem];
    DataType* data = new DataType[num_elem];

    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ElementCentroids(num_pos,pos,num_elem,conn,centroids);

    int Time = CurrentTime(base::GH(),0);
    double t = GetPhysicalTime(IOfunc,Time,0);
    Interpolation().PointsValues(IOfunc,t,num_elem,reinterpret_cast<const point_type*>(centroids),data);
    Interpolation().ArrayCombine(VizServer,num_elem,data);

    geom_point_type sum(0.);
    if (MY_PROC==VizServer) {
      DataType* area = new DataType[num_elem];
      geom_point_type* normal = new geom_point_type[num_elem];

      ops::ElementAreas(num_pos,pos,num_elem,conn,area);
      ops::ElementNormals(num_pos,pos,num_elem,conn,normal);     
      for (register int i=0; i<num_elem; i++) 
        if (data[i] != Interpolation().ErrorValue())
          sum += data[i]*area[i]*normal[i];

      delete [] area;
      delete [] normal;
    }

    delete [] centroids;
    delete [] data;

    return sum;
  } 

  geom_point_type IntegrateVector(const geom_point_type* data, 
                                  const int num_pos, const DataType* vertices,
                                  const int num_elem, const int* connections) {
    geom_point_type sum(0.);
    DataType* area = new DataType[num_elem];
    geom_point_type* normal = new geom_point_type[num_elem];
    
    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ops::ElementAreas(num_pos,pos,num_elem,conn,area);   
    for (register int i=0; i<num_elem; i++) 
      if (data[i][0] != Interpolation().ErrorValue())
        sum += area[i]*data[i];
    
    delete [] area;

    return sum;
  } 

  DataType IntegrateVectorNormal(const geom_point_type* data, 
                                 const int num_pos, const DataType* vertices,
                                 const int num_elem, const int* connections) {
    DataType sum(0.);
    DataType* area = new DataType[num_elem];
    geom_point_type* normal = new geom_point_type[num_elem];
    
    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ops::ElementAreas(num_pos,pos,num_elem,conn,area);
    ops::ElementNormals(num_pos,pos,num_elem,conn,normal);  
    for (register int i=0; i<num_elem; i++) 
      if (data[i][0] != Interpolation().ErrorValue())
        sum += area[i]*geom::computeDotProduct(data[i],normal[i]);
    
    delete [] area;
    delete [] normal;
    
    return sum;
  } 

  geom_point_type IntegrateTensorNormal(const geom_point_type* data, 
                                        const int num_pos, const DataType* vertices,
                                        const int num_elem, const int* connections) {
    geom_point_type sum(0.);
    DataType* area = new DataType[num_elem];
    geom_point_type* normal = new geom_point_type[num_elem];
    
    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ops::ElementAreas(num_pos,pos,num_elem,conn,area);
    ops::ElementNormals(num_pos,pos,num_elem,conn,normal);
    
    for (register int i=0; i<num_elem; i++) 
      for (register int d=0; d<dim; d++) 
      if (data[dim*i+d][0] != Interpolation().ErrorValue())
        sum[d] += area[i]*geom::computeDotProduct(data[dim*i+d],normal[i]);
    
    delete [] area;
    delete [] normal;
    
    return sum;
  } 

  geom_point_type IntegrateMoment(const DataType* data, 
                                  const int num_pos, const DataType* vertices,
                                  const int num_elem, const int* connections,
                                  const geom_point_type& point) {
    geom_point_type sum(0.);
    DataType* area = new DataType[num_elem];
    geom_point_type* normal = new geom_point_type[num_elem];
    geom_point_type* cross = new geom_point_type[num_elem];
    geom_point_type* centroids = new geom_point_type[num_elem];
    
    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ElementCentroids(num_pos,pos,num_elem,conn,centroids);
    ops::ElementAreas(num_pos,pos,num_elem,conn,area);
    ops::ElementNormals(num_pos,pos,num_elem,conn,normal);
    
    for (register int i=0; i<num_elem; i++)
      if (data[i] != Interpolation().ErrorValue()) {
        centroids[i] -= point;
        normal[i] *= data[i];
      }
    ops::CrossProducts(num_elem,centroids,normal,cross);
    for (register int i=0; i<num_elem; i++) 
      if (data[i] != Interpolation().ErrorValue()) 
        sum += area[i]*cross[i];
    
    delete [] area;
    delete [] normal;
    delete [] cross;
    delete [] centroids;

    return sum;
  }

  geom_point_type IntegrateMoment(grid_fct_type& IOfunc, 
                                  const int num_pos, const DataType* vertices,
                                  const int num_elem, const int* connections,
                                  const geom_point_type& point) {
    geom_point_type* centroids = new geom_point_type[num_elem];
    DataType* data = new DataType[num_elem];

    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ElementCentroids(num_pos,pos,num_elem,conn,centroids);

    int Time = CurrentTime(base::GH(),0);
    double t = GetPhysicalTime(IOfunc,Time,0);
    Interpolation().PointsValues(IOfunc,t,num_elem,reinterpret_cast<const point_type*>(centroids),data);
    Interpolation().ArrayCombine(VizServer,num_elem,data);
 
    geom_point_type sum(0.);
    if (MY_PROC==VizServer) {
      DataType* area = new DataType[num_elem];
      geom_point_type* normal = new geom_point_type[num_elem];
      geom_point_type* cross = new geom_point_type[num_elem];

      ops::ElementAreas(num_pos,pos,num_elem,conn,area);
      ops::ElementNormals(num_pos,pos,num_elem,conn,normal);
      for (register int i=0; i<num_elem; i++)
        if (data[i] != Interpolation().ErrorValue()) {
          centroids[i] -= point;
          normal[i] *= data[i];
        }
      ops::CrossProducts(num_elem,centroids,normal,cross);
      for (register int i=0; i<num_elem; i++) 
        if (data[i] != Interpolation().ErrorValue()) 
          sum += area[i]*cross[i];

      delete [] area;
      delete [] normal;
      delete [] cross;
    }

    delete [] centroids;
    delete [] data;

    return sum;
  }

  geom_point_type IntegrateVectorMoment(const geom_point_type* data, 
                                        const int num_pos, const DataType* vertices,
                                        const int num_elem, const int* connections,
                                        const geom_point_type& point) {
    geom_point_type sum(0.);
    DataType* area = new DataType[num_elem];
    geom_point_type* cross = new geom_point_type[num_elem];
    geom_point_type* centroids = new geom_point_type[num_elem];
    
    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ElementCentroids(num_pos,pos,num_elem,conn,centroids);
    ops::ElementAreas(num_pos,pos,num_elem,conn,area);    
    for (register int i=0; i<num_elem; i++)
      if (data[i][0] != Interpolation().ErrorValue()) 
        centroids[i] -= point;
    ops::CrossProducts(num_elem,centroids,data,cross);
    for (register int i=0; i<num_elem; i++) 
      if (data[i][0] != Interpolation().ErrorValue()) 
        sum += area[i]*cross[i];
    
    delete [] area;
    delete [] cross;
    delete [] centroids;

    return sum;
  }

  geom_point_type IntegrateTensorMoment(const geom_point_type* data, 
                                        const int num_pos, const DataType* vertices,
                                        const int num_elem, const int* connections,
                                        const geom_point_type& point) {
    geom_point_type sum(0.);
    DataType* area = new DataType[num_elem];
    geom_point_type* normal = new geom_point_type[num_elem];
    geom_point_type* cross = new geom_point_type[num_elem];
    geom_point_type* datavec = new geom_point_type[num_elem];
    geom_point_type* centroids = new geom_point_type[num_elem];
    
    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ElementCentroids(num_pos,pos,num_elem,conn,centroids);
    ops::ElementAreas(num_pos,pos,num_elem,conn,area);
    ops::ElementNormals(num_pos,pos,num_elem,conn,normal);
    for (register int i=0; i<num_elem; i++) {
      centroids[i] -= point;
      for (register int d=0; d<dim; d++) 
        if (data[dim*i+d][0] != Interpolation().ErrorValue())
          datavec[i][d] = geom::computeDotProduct(data[dim*i+d],normal[i]);
        else
          datavec[i][d] = Interpolation().ErrorValue();
    }
    ops::CrossProducts(num_elem,centroids,datavec,cross);
    for (register int i=0; i<num_elem; i++) 
      if (datavec[i][0] != Interpolation().ErrorValue()) 
        sum += area[i]*cross[i];
    
    delete [] area;
    delete [] normal;
    delete [] cross;
    delete [] datavec;
    delete [] centroids;

    return sum;
  }

  void ElementCentroids(const int num_pos, const geom_point_type* pos,
                        const int num_elem, const multi_index_type* conn,
                        geom_point_type* centroids) {
    for (register int i=0; i<num_elem; i++) {
      const multi_index_type& face = conn[i];
      centroids[i] = pos[face[0]];
      for (register int d=1; d<dim; d++)
        centroids[i] += pos[face[d]];
      centroids[i] /= dim;
    }
  }

  void ElementCentroids(const int num_pos, const DataType* vertices,
                        const int num_elem, const int* connections,
                        geom_point_type* centroids) {
    const geom_point_type* pos = reinterpret_cast<const geom_point_type*>(vertices);
    const multi_index_type* conn = reinterpret_cast<const multi_index_type*>(connections);
    ElementCentroids(num_pos,pos,num_elem,conn,centroids);
  }
  
  inline void SetNcnt(const int& cnt) { _Ncnt = cnt; }
  inline const int& Ncnt() const { return _Ncnt; }
  inline const int& OutputType() const { return _OutputType; }
  inline const std::string& OutputName(const int cnt) const 
  { assert (cnt>=0 && cnt<MAXCOMPONENTS); return CompName[cnt]; }

  inline void SetInterpolation(interpolation_type* interpol) { _Interpolation = interpol; }
  inline interpolation_type& Interpolation() { return *_Interpolation; }
  inline const interpolation_type& Interpolation() const { return *_Interpolation; }

protected:
  interpolation_type* _Interpolation;
  int _Ncnt, _OutputType;
  int StepFormat, StepDirectories;
  std::string BaseFileName;
  std::string* CompName;
  std::string CompDir;
  ControlDevice LocCtrl;
};


#endif

---