Maxwell-TD-Scattering/MaxwellTD_GPU/fem/plane_wave_mesh.cc

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <string>
#include <iomanip>

#include "plane_wave_mesh.h"
#include "femgrp.h"
#include "face.h"
#include "node.h"

PlaneWaveMesh::PlaneWaveMesh(){
  memset(PW_Name, 0, 180 * sizeof(char));
  theta_ = 0.;
  phi_ = 0.;
  imp_ = 0.;
  faceCNT = 0;
  objCNT = 0;
  fcArray = 0;
  ndArray = 0;
  objMAP = 0;
}

PlaneWaveMesh::~PlaneWaveMesh(){
  delete [] ndArray;
  delete [] fcArray;
  delete [] objMAP;
}

void PlaneWaveMesh::setFaceCnt(int faceCnt){
  faceCNT = faceCnt;
  fcArray = new face*[faceCNT];
  memset(fcArray, 0, faceCNT * sizeof(face*));
}

void PlaneWaveMesh::setNodeCnt(int nodeCnt){
  nodeCNT = nodeCnt;
  ndArray = new node*[nodeCNT];
  memset(ndArray, 0, nodeCNT * sizeof(node*));
}

int PlaneWaveMesh::getObjMap(int obj){
  int i;
  for(i = 0; i < objCNT; i++){
    if(objMAP[i] == obj) 
      return i;
  }
  return -1;
}

void PlaneWaveMesh::insertObj(int obj){
  int i;
  for(i = 0; i < objCNT; i++){
    if(objMAP[i] == obj) 
      return;
  }
  objMAP[objCNT] = obj;
  objCNT++;
}

void PlaneWaveMesh::writeVtk(char* fileName){
  // int i, j;
  // char hFileName[STRLEN];
  // char eFileName[STRLEN];
  // FILE* fdE;
  // FILE* fdH;
  // char cc;
  // int dummy;
  // vtr dummyVtr;

  // // fill the port field with averaged values
  // vtr* eField = new vtr[nodeCNT];
  // vtr* hField = new vtr[nodeCNT];
  // int*  count = new int[nodeCNT];
  // memset(count, 0, nodeCNT*sizeof(int));

  // int index;
  // cVtr evtr;
  // cVtr hvtr;
  // vtr  tmpE; // change here
  // vtr  tmpH;
  // for(i = 0; i < faceCNT; i++){
  //   for(j = 0; j < SEC_ORDER_FACE_NODES; j++){
  //     fread(&evtr, sizeof(cVtr), 1, fdE);
  //     fread(&hvtr, sizeof(cVtr), 1, fdH);
  //     if(j < 3){
  //       index = globToLocMap_->find(fcArray[i]->getNode(j)->getid())->second;
  //       tmpE.setvtr(evtr.x.real(), evtr.y.real(), evtr.z.real());
  //       tmpH.setvtr(hvtr.x.real(), hvtr.y.real(), hvtr.z.real());
  //       eField[index] = eField[index] + tmpE;
  //       hField[index] = hField[index] + tmpH;
  //       (count[index])++;
  //     }
  //   }
  // }

  // for(i = 0; i < nodeCNT; i++){
  //   eField[i] = eField[i] / static_cast<fp_t>(count[i]);
  //   hField[i] = hField[i] / static_cast<fp_t>(count[i]);
  // }

  // node** locNodeArray = new node*[nodeCNT];
  // for(i = 0; i < nodeCNT; i++){
  //   node& Node = *(ndArray[i]);
  //   int index = globToLocMap_->find(Node.getid())->second;
  //   locNodeArray[index] = new node(index, Node.getPType(), Node.getCoord().getx(), Node.getCoord().gety(), Node.getCoord().getz());
  // }

  // face** locFaceArray = new face*[faceCNT];
  // for(i = 0; i < faceCNT; i++){
  //   face& Face = *(fcArray[i]);
  //   locFaceArray[i] = new face(Face);
  //   locFaceArray[i]->setFace(
  //       locNodeArray[globToLocMap_->find(Face.getNode(0)->getid())->second],
  //       locNodeArray[globToLocMap_->find(Face.getNode(1)->getid())->second],
  //       locNodeArray[globToLocMap_->find(Face.getNode(2)->getid())->second]);
  // }

  // VtkWriter vtkWriter(1.);
  // vtkWriter.writeTriUg(fileName, nodeCNT, locNodeArray, faceCNT, locFaceArray, eField, hField, 1);

  // for(i = 0; i < nodeCNT; i++) 
  //   delete locNodeArray[i];
  // delete [] locNodeArray;
  // for(i = 0; i < faceCNT; i++) 
  //   delete locFaceArray[i];
  // delete [] locFaceArray;
  // delete [] eField;
  // delete [] hField;
  // delete [] count;
}


void PlaneWaveMesh::computeBoundingBox() 
{
  if (nodeCNT <= 0 || !ndArray) {
    std::cerr << "[PlaneWaveMesh::computeBoundingBox] No nodes available.\n";
    return;
  }

  // Initialize with the first node (non-parallel)
  const vtr& coord = ndArray[0]->getCoord();
  xmin = xmax = coord.getx();
  ymin = ymax = coord.gety();
  zmin = zmax = coord.getz();

  // Parallel reduction for bounding box
  // The reduction(min: ...) and reduction(max: ...) clauses are used to ensure that min/max values are correctly computed in parallel.
  // Each thread should compute a local copy of this variable, and at the end of the loop, combine all local results into one.
  #pragma omp parallel for reduction(min:xmin, ymin, zmin) reduction(max:xmax, ymax, zmax)
  for (int i = 1; i < nodeCNT; ++i) {
    const vtr& coord = ndArray[i]->getCoord();
    const fp_t x = coord.getx();
    const fp_t y = coord.gety();
    const fp_t z = coord.getz();

    if (x < xmin) xmin = x;
    if (x > xmax) xmax = x;

    if (y < ymin) ymin = y;
    if (y > ymax) ymax = y;

    if (z < zmin) zmin = z;
    if (z > zmax) zmax = z;
  }

  std::cout << "Bounding Box for Plane Wave Mesh Nodes:\n";
  std::cout << "  x: [" << xmin << ", " << xmax << "]\n";
  std::cout << "  y: [" << ymin << ", " << ymax << "]\n";
  std::cout << "  z: [" << zmin << ", " << zmax << "]\n";
}