neBEMInterface.h File Reference

#include <stdio.h>
#include <time.h>
#include "Vector.h"

Go to the source code of this file.

Macros
#define	INTFACEGLOBAL   extern
#define	neBEMtwopi   6.283185307179586476925287
#define	neBEMrtod   57.2957795
#define	MODULUS(p)

Functions
INTFACEGLOBAL double	neBEMTimeElapsed (clock_t, clock_t)
INTFACEGLOBAL int	neBEMMessage (const char *message)
INTFACEGLOBAL int	neBEMGetNbOfLines (const char *fname)
INTFACEGLOBAL double	neBEMChkInPoly (int nvert, Point3D node[], Point3D pt)
INTFACEGLOBAL int	neBEMSetDefaults (void)
INTFACEGLOBAL int	ReadInitFile (char filename[])
INTFACEGLOBAL int	neBEMGetInputsFromFiles (void)
INTFACEGLOBAL int	neBEMInitialize (void)
INTFACEGLOBAL int	neBEMReadGeometry (void)
INTFACEGLOBAL int	neBEMGetNbPrimitives (void)
INTFACEGLOBAL int	neBEMGetPrimitive (int prim, int *nvertex, double xvert[], double yvert[], double zvert[], double *xnorm, double *ynorm, double *znorm, int *volref1, int *volref2)
INTFACEGLOBAL int	neBEMVolumeDescription (int ivol, int *shape, int *material, double *epsilon, double *potential, double *charge, int *boundarytype)
INTFACEGLOBAL int	neBEMVolumePoint (double x, double y, double z)
INTFACEGLOBAL void	neBEMVolumePrimitives (int volume, int *nprim, int primlist[])
INTFACEGLOBAL int	neBEMGetPeriodicities (int prim, int *ix, int *jx, double *sx, int *iy, int *jy, double *sy, int *iz, int *jz, double *sz)
INTFACEGLOBAL int	neBEMGetMirror (int prim, int *ix, int *jx, double *sx, int *iy, int *jy, double *sy, int *iz, int *jz, double *sz)
INTFACEGLOBAL int	neBEMGetBoundingPlanes (int *ixmin, double *cxmin, double *vxmin, int *ixmax, double *cxmax, double *vxmax, int *iymin, double *cymin, double *vymin, int *iymax, double *cymax, double *vymax, int *izmin, double *czmin, double *vzmin, int *izmax, double *czmax, double *vzmax)
INTFACEGLOBAL int	neBEMDiscretize (int **elementNbs)
INTFACEGLOBAL int	neBEMBoundaryInitialConditions (void)
INTFACEGLOBAL int	neBEMSolve (void)
INTFACEGLOBAL int	neBEMPF (Point3D *, double *, Vector3D *)
INTFACEGLOBAL double	neBEMVolumeCharge (int volume)
INTFACEGLOBAL int	neBEMPrepareWeightingField (int NbPrimsWtField, int PrimListWtField[])
INTFACEGLOBAL void	neBEMDeleteWeightingField (int IdWtField)
INTFACEGLOBAL void	neBEMDeleteAllWeightingFields (void)
INTFACEGLOBAL double	neBEMWeightingField (Point3D *point, Vector3D *field, int IdWtField)
INTFACEGLOBAL int	neBEMEnd (void)
INTFACEGLOBAL int	CreateDirStr (void)
INTFACEGLOBAL int	WritePrimitives (void)
INTFACEGLOBAL int	ReadPrimitives (void)
INTFACEGLOBAL int	WriteElements (void)
INTFACEGLOBAL int	ReadElements (void)

Variables
INTFACEGLOBAL clock_t	startClock
INTFACEGLOBAL clock_t	stopClock
INTFACEGLOBAL int	neBEMState
INTFACEGLOBAL int	OptDeviceFile
INTFACEGLOBAL char	DeviceInputFile [256]
INTFACEGLOBAL int	NbThreads
INTFACEGLOBAL int	OptPrintPrimaryDetails
INTFACEGLOBAL int	OptPrintVolumeDetails
INTFACEGLOBAL int	OptGnuplot
INTFACEGLOBAL int	OptGnuplotPrimitives
INTFACEGLOBAL int	OptGnuplotElements
INTFACEGLOBAL FILE *	fgnuPrim
INTFACEGLOBAL FILE *	fgnuElem
INTFACEGLOBAL FILE *	fgnuMesh
INTFACEGLOBAL int	OptPrintVertexAndNormal
INTFACEGLOBAL int	OptPrimitiveFiles
INTFACEGLOBAL int	OptElementFiles
INTFACEGLOBAL int	neBEMPFCallCntr
INTFACEGLOBAL int	OptReuseDir

Macro Definition Documentation

INTFACEGLOBAL

#define INTFACEGLOBAL   extern

Definition at line 7 of file neBEMInterface.h.

Referenced by neBEMChkInPoly(), neBEMGetBoundingPlanes(), neBEMGetInputsFromFiles(), neBEMGetMirror(), neBEMGetNbOfLines(), neBEMGetNbPrimitives(), neBEMGetPeriodicities(), neBEMGetPrimitive(), neBEMSetDefaults(), neBEMVolumeDescription(), neBEMVolumePoint(), neBEMVolumePrimitives(), and ReadInitFile().

MODULUS

#define MODULUS

(

p

)

Value:

(sqrt(p.X * p.X + p.Y * p.Y + p.Z * p.Z))

Definition at line 35 of file neBEMInterface.h.

Referenced by neBEMChkInPoly().

neBEMrtod

#define neBEMrtod   57.2957795

Definition at line 34 of file neBEMInterface.h.

neBEMtwopi

#define neBEMtwopi   6.283185307179586476925287

Definition at line 33 of file neBEMInterface.h.

Referenced by InitChargingUp(), and neBEMChkInPoly().

Function Documentation

CreateDirStr()

INTFACEGLOBAL int CreateDirStr

(

void

)

Definition at line 2702 of file neBEMInterface.c.

                       {
  char strModelCntr[10], strMeshCntr[10], strBCCntr[10], strPPCntr[10];
  int CreateOrUseDir(char[]);
  int CreateDirOrQuit(char[]);
 
  snprintf(strModelCntr, 10, "/Model%d", ModelCntr);
  snprintf(strMeshCntr, 10, "/M%d", MeshCntr);
  snprintf(strBCCntr, 10, "/BC%d", BCCntr);
  snprintf(strPPCntr, 10, "/PP%d", PPCntr);
 
  strcpy(ModelOutDir, DeviceOutDir);
  strcat(ModelOutDir, strModelCntr);
  strcpy(NativeOutDir, ModelOutDir);
  strcat(NativeOutDir, "/neBEMNatives/");
  strcpy(NativePrimDir, NativeOutDir);
  strcat(NativePrimDir, "Primitives/");
  strcpy(MeshOutDir, ModelOutDir);
  strcat(MeshOutDir, strMeshCntr);
  strcpy(BCOutDir, MeshOutDir);
  strcat(BCOutDir, strBCCntr);
  strcpy(PPOutDir, BCOutDir);
  strcat(PPOutDir, strPPCntr);
 
  // Create DeviceOutDir, if necessary
  int fstatus = CreateOrUseDir(DeviceOutDir);
  if (fstatus != 0) {
    neBEMMessage("CreateDirStr - CreateOrUseDir");
    return -1;
  }
 
  if (NewModel) {
    // create ModelOutDir
    if (OptReuseDir) {
      fstatus = CreateOrUseDir(ModelOutDir);
      fstatus = CreateOrUseDir(NativeOutDir);
      fstatus = CreateOrUseDir(NativePrimDir);
    } else {
      fstatus = CreateDirOrQuit(ModelOutDir);
      fstatus = CreateDirOrQuit(NativeOutDir);
      fstatus = CreateDirOrQuit(NativePrimDir);
    }
    if (fstatus != 0) {
      neBEMMessage("CreateDirStr - ModelOutDir");
      return -1;
    }
  }
 
  if (NewMesh) {
    // create MeshOutDir
    if (OptReuseDir)
      fstatus = CreateOrUseDir(MeshOutDir);
    else
      fstatus = CreateDirOrQuit(MeshOutDir);
    if (fstatus != 0) {
      neBEMMessage("CreateDirStr - MeshOutDir");
      return -1;
    }
  }
 
  if (NewBC) {
    // create BCOutDir
    if (OptReuseDir)
      fstatus = CreateOrUseDir(BCOutDir);
    else
      fstatus = CreateDirOrQuit(BCOutDir);
    if (fstatus != 0) {
      neBEMMessage("CreateDirStr - BCOutDir");
      return -1;
    }
  }
 
  if (NewPP) {
    // create PPOutDir
    if (OptReuseDir)
      fstatus = CreateOrUseDir(PPOutDir);
    else
      fstatus = CreateDirOrQuit(PPOutDir);
    if (fstatus != 0) {
      neBEMMessage("CreateDirStr - PPOutDir");
      return -1;
    }
  }
 
  // Create other relevant sub-directories
  char subdir[256];
 
  strcpy(subdir, ModelOutDir);
  strcat(subdir, "/Primitives/");
  if (OptReuseDir)
    fstatus = CreateOrUseDir(subdir);
  else
    fstatus = CreateDirOrQuit(subdir);
 
  strcpy(subdir, MeshOutDir);
  strcat(subdir, "/Elements/");
  if (OptReuseDir)
    fstatus = CreateOrUseDir(subdir);
  else
    fstatus = CreateDirOrQuit(subdir);
 
  strcpy(subdir, MeshOutDir);
  strcat(subdir, "/GViewDir/");
  if (OptReuseDir)
    fstatus = CreateOrUseDir(subdir);
  else
    fstatus = CreateDirOrQuit(subdir);
 
  return (0);
}  // CreateDirStr ends

Referenced by neBEMInitialize().

neBEMBoundaryInitialConditions()

INTFACEGLOBAL int neBEMBoundaryInitialConditions

(

void

)

Definition at line 1756 of file neBEMInterface.c.

                                         {
  startClock = clock();
 
  // The boundary conditions can be set only with neBEMState == 4 or 7
  // This state is assigned either after element discretization has been
  // completed or in a condition when we are looking for modifying only the
  // boundary condition for a device having same geometry (hence, the same
  // inverted influence coefficient matrix)
  if ((neBEMState == 4) || (neBEMState == 7)) {
    int fstatus = BoundaryConditions();
    if (fstatus != 0) {
      neBEMMessage("neBEMBondaryInitialConditions - BoundaryConditions");
      return -1;
    }
    fstatus = InitialConditions();
    if (fstatus != 0) {
      neBEMMessage("neBEMBondaryInitialConditions - InitialConditions");
      return -1;
    }
    if (neBEMState == 4) neBEMState = 5;  // create LHMatrix, invert etc
    if (neBEMState == 7) neBEMState = 8;  // assume LHS and inversion to be over
  } else {
    printf("Boundary & initial conditions can be set in states 4 / 7 ...\n");
    printf("returning ...\n");
    return (-1);
  }
 
  stopClock = clock();
  neBEMTimeElapsed(startClock, stopClock);
  printf("to setup boundary and initial conditions.\n");
 
  return (0);
}  // neBEMBoundaryInitialConditions ends

neBEMChkInPoly()

INTFACEGLOBAL double neBEMChkInPoly	(	int	nvert,
		Point3D	node[],
		Point3D	pt )

neBEMDeleteAllWeightingFields()

INTFACEGLOBAL void neBEMDeleteAllWeightingFields

(

void

)

Definition at line 2593 of file neBEMInterface.c.

                                         {
  const int MaxWtField = MAXWtFld;  // being used while allocating memory
  for (int id = 1; id < MaxWtField; ++id) {  // count from 1
    free(WtFieldChDen[id]);
    free(AvWtChDen[id]);
  }
  free(WtFieldChDen);
  free(AvWtChDen);
}

neBEMDeleteWeightingField()

INTFACEGLOBAL void neBEMDeleteWeightingField

(

int

IdWtField

)

Definition at line 2587 of file neBEMInterface.c.

                                              {
  free(WtFieldChDen[IdWtField]);
  free(AvWtChDen[IdWtField]);
}

neBEMDiscretize()

INTFACEGLOBAL int neBEMDiscretize

(

int **

elementNbs

)

Definition at line 1454 of file neBEMInterface.c.

                                               {
  // int fstatus;
 
  // For a model and a mesh that were defined before and for which data were
  // stored in two files - one for primitives and one for elements
  // The following operation does not assume that the primitives have been read
  // in from a stored file. In essence, these two read-in-s are maintained
  // independent of each other. However, for a stored element file to be useful,
  // both the model and mesh have to be old (not new).
  if ((!NewModel) && (!NewMesh) && (!NewBC) && (OptStoreElements)) {
    int fstatus = ReadElements();
    if (fstatus) {
      neBEMMessage("neBEMDiscretize - problem reading stored Elements.\n");
      return -1;
    }
    neBEMState = 4;
    return 0;
  }
 
  // Otherwise, continue with fresh discretization
  if (neBEMState != 3) {
    printf("discretization can continue only in State 3 ...\n");
    return -1;
  }
 
  // Only the number of primitives has been ascertained.
  // All the rest globally important numbers will be determined in this function
  NbSurfs = 0;
  NbWires = 0;
  NbElements = 0;
 
  // Here, the primitive can be analyzed and the elements necessary to
  // discretize it, may be determined. A user hint may help that can be supplied
  // during the setting up of the device. The hint may be as naive as gross,
  // coarse, medium, regular, fine depending on which the element size (in %
  // of the device / primitive) may be decided upon.
  // Earlier method of specifying the number of primitives on a priori has been
  // over-ridden.
  // Now we prescribe the length / area of each element and discretize each
  // primitive such that each element has an length / area close to the
  // suggested value.
  // Since the number of elements are being decided here, all the shapes will
  // have to be one among wire, right triangle and restangle. Any decomposition
  // of arbitrary polygons into rectangles (as many as possible) and right
  // triangles will have to be done earlier. All the counts have been
  // incremented by one to be on the safe side! The additional memory can be
  // minimized through a more careful computation of the required allocation for
  // each type of primitive.
  char MeshLogFile[256];
 
  strcpy(MeshLogFile, MeshOutDir);
  strcat(MeshLogFile, "/MeshLog.out");
  fMeshLog = fopen(MeshLogFile, "w");
  fprintf(fMeshLog, "Details of primitive discretization\n");
 
  for (int prim = 1; prim <= NbPrimitives; ++prim) {
    if (NbVertices[prim] == 4) {
      NbSurfSegX[prim] = NbElemsOnPrimitives[prim][1];
      NbSurfSegZ[prim] = NbElemsOnPrimitives[prim][2];
      int fstatus =
          AnalyzePrimitive(prim, &NbSurfSegX[prim], &NbSurfSegZ[prim]);
      if (fstatus == 0) {
        neBEMMessage("neBEMDiscretize - AnalyzePrimitve");
        return -1;
      }
      NbElements += (NbSurfSegX[prim] + 1) * (NbSurfSegZ[prim] + 1);
    }
    if (NbVertices[prim] == 3) {
      NbSurfSegX[prim] = NbElemsOnPrimitives[prim][1];
      NbSurfSegZ[prim] = NbElemsOnPrimitives[prim][2];
      int fstatus =
          AnalyzePrimitive(prim, &NbSurfSegX[prim], &NbSurfSegZ[prim]);
      if (fstatus == 0) {
        neBEMMessage("neBEMDiscretize - AnalyzePrimitive");
        return -1;
      }
      NbElements += (NbSurfSegX[prim] + 1) * (NbSurfSegZ[prim] + 1);
    }
    if (NbVertices[prim] == 2) {
      int itmp;
      NbWireSeg[prim] = NbElemsOnPrimitives[prim][1];
      int fstatus = AnalyzePrimitive(prim, &NbWireSeg[prim], &itmp);
      if (fstatus == 0) {
        neBEMMessage("neBEMDiscretize - AnalyzePrimitive");
        return -1;
      }
      NbElements += (NbWireSeg[prim] + 1);
    }
    if (DebugLevel == 101) {
      if (NbVertices[prim] == 2) {
        printf("Primitive %d to be discretized into %d elements.\n", prim,
               NbWireSeg[prim]);
      } else {
        printf("Primitive %d to be discretized into %d X %d elements.\n", prim,
               NbSurfSegX[prim], NbSurfSegZ[prim]);
      }
    }
  }
  printf("Memory allocated for maximum %d elements.\n", NbElements);
  fclose(fMeshLog);
 
  // Allocate enough space to store all the elements
  if (neBEMState == 3) {
    printf("neBEMDiscretize: NbElements = %d, sizeof(Element) = %zu\n",
           NbElements, sizeof(Element));
    if (EleArr) {
      Element *tmp = (Element *)realloc(EleArr, NbElements * sizeof(Element));
      if (tmp != NULL) {
        EleArr = tmp;
        EleCntr = 0;
      } else {
        free(EleArr);
        printf("neBEMDiscretize: Re-allocating EleArr failed.\n");
        return (1);
      }
      printf("neBEMDiscretize: Re-allocated EleArr.\n");
    }  // if EleArr => re-allocation
    else {
      EleArr = (Element *)malloc(NbElements * sizeof(Element));
      if (EleArr == NULL) {
        neBEMMessage("neBEMDiscretize - EleArr malloc");
        return -1;
      }
    }  // else EleArr => fresh allocation
  }    // neBEMState == 3
 
  // Prepare a data file that will contain the plotting information of the
  // the primitives and the elements
  if (OptGnuplot) {
    char GnuFile[256];
    strcpy(GnuFile, MeshOutDir);
    strcat(GnuFile, "/GViewDir/gPrimView.gp");
    fgnuPrim = fopen(GnuFile, "w");
    fprintf(fgnuPrim, "set title \"neBEM primitives in gnuplot VIEWER\"\n");
    // fprintf(fgnu, "#set label 1 \'LengthScale = %d\', LengthScale, right\n");
    fprintf(fgnuPrim, "#set pm3d\n");
    fprintf(fgnuPrim, "#set style data pm3d\n");
    fprintf(fgnuPrim, "#set palette model CMY\n");
    fprintf(fgnuPrim, "set hidden3d\n");
    fprintf(fgnuPrim, "set nokey\n");
    fprintf(fgnuPrim, "set xlabel \"X\"\n");
    fprintf(fgnuPrim, "set ylabel \"Y\"\n");
    fprintf(fgnuPrim, "set zlabel \"Z\"\n");
    fprintf(fgnuPrim, "set view 70, 335, 1, 1\n");
    fprintf(fgnuPrim, "\nsplot \\\n");
 
    strcpy(GnuFile, MeshOutDir);
    strcat(GnuFile, "/GViewDir/gElemView.gp");
    fgnuElem = fopen(GnuFile, "w");
    fprintf(fgnuElem, "set title \"neBEM elements in gnuplot VIEWER\"\n");
    // fprintf(fgnu, "#set label 1 \'LengthScale = %d\', LengthScale, right\n");
    fprintf(fgnuElem, "#set pm3d\n");
    fprintf(fgnuElem, "#set style data pm3d\n");
    fprintf(fgnuElem, "#set palette model CMY\n");
    fprintf(fgnuElem, "set hidden3d\n");
    fprintf(fgnuElem, "set nokey\n");
    fprintf(fgnuElem, "set xlabel \"X\"\n");
    fprintf(fgnuElem, "set ylabel \"Y\"\n");
    fprintf(fgnuElem, "set zlabel \"Z\"\n");
    fprintf(fgnuElem, "set view 70, 335, 1, 1\n");
    fprintf(fgnuElem, "\nsplot \\\n");
 
    strcpy(GnuFile, MeshOutDir);
    strcat(GnuFile, "/GViewDir/gMeshView.gp");
    fgnuMesh = fopen(GnuFile, "w");
    fprintf(fgnuMesh, "set title \"neBEM mesh in gnuplot VIEWER\"\n");
    // fprintf(fgnu, "#set label 1 \'LengthScale = %d\', LengthScale, right\n");
    fprintf(fgnuMesh, "#set pm3d\n");
    fprintf(fgnuMesh, "#set style data pm3d\n");
    fprintf(fgnuMesh, "#set palette model CMY\n");
    fprintf(fgnuMesh, "set hidden3d\n");
    fprintf(fgnuMesh, "set nokey\n");
    fprintf(fgnuMesh, "set xlabel \"X\"\n");
    fprintf(fgnuMesh, "set ylabel \"Y\"\n");
    fprintf(fgnuMesh, "set zlabel \"Z\"\n");
    fprintf(fgnuMesh, "set view 70, 335, 1, 1\n");
    fprintf(fgnuMesh, "\nsplot \\\n");
  }
 
  for (int prim = 1; prim <= NbPrimitives; ++prim) {
    switch (PrimType[prim]) {
      int fstatus;
      case 3:  // triangular surface
      case 4:  // rectangular surface
        ++NbSurfs;
        fstatus = SurfaceElements(
            prim, NbVertices[prim], XVertex[prim], YVertex[prim], ZVertex[prim],
            XNorm[prim], YNorm[prim], ZNorm[prim], VolRef1[prim], VolRef2[prim],
            InterfaceType[prim], ApplPot[prim], ApplCh[prim], Lambda[prim],
            NbSurfSegX[prim], NbSurfSegZ[prim]);
        if (fstatus != 0) {
          neBEMMessage("neBEMDiscretize - SurfaceElements");
          return -1;
        }
        break;
      case 2:       // wire - a wire presumably has only 2 vertices
        ++NbWires;  // it has one radius and one segmentation information
        fstatus = WireElements(
            prim, NbVertices[prim], XVertex[prim], YVertex[prim], ZVertex[prim],
            Radius[prim], VolRef1[prim], VolRef2[prim], InterfaceType[prim],
            ApplPot[prim], ApplCh[prim], Lambda[prim], NbWireSeg[prim]);
        if (fstatus != 0) {
          neBEMMessage("neBEMDiscretize - WireElements");
          return -1;
        }
        break;
 
      default:
        printf("PrimType out of range in CreateElements ... exiting ...\n");
        exit(-1);
    }  // switch PrimType ends
  }    // loop on prim number ends
 
  if (OptGnuplot) {
    fprintf(fgnuPrim, "\n\npause-1");
    fclose(fgnuPrim);
    fprintf(fgnuElem, "\n\npause-1");
    fclose(fgnuElem);
    fprintf(fgnuMesh, "\n\npause-1");
    fclose(fgnuMesh);
  }
 
  // If the required memory exceeds the maximum allowed number of elements
  if (EleCntr > NbElements) {
    neBEMMessage("neBEMDiscretize - EleCntr more than NbElements!");
    return -1;
  }
 
  // Check whether collocation points overlap
  {
    int startcntr = 1, cntr1, cntr2;
    Point3D pt1, pt2;
    double dist;
    for (cntr1 = startcntr; cntr1 <= EleCntr; ++cntr1) {
      pt1.X = (EleArr + cntr1 - 1)->BC.CollPt.X;
      pt1.Y = (EleArr + cntr1 - 1)->BC.CollPt.Y;
      pt1.Z = (EleArr + cntr1 - 1)->BC.CollPt.Z;
      for (cntr2 = cntr1 + 1; cntr2 <= EleCntr; ++cntr2) {
        pt2.X = (EleArr + cntr2 - 1)->BC.CollPt.X;
        pt2.Y = (EleArr + cntr2 - 1)->BC.CollPt.Y;
        pt2.Z = (EleArr + cntr2 - 1)->BC.CollPt.Z;
 
        dist = GetDistancePoint3D(&pt1, &pt2);
        if (dist <= MINDIST)  {
          // we need a linked-list here so that the overlapped
          // element is easily deleted and the rest upgraded immediately
 
          // Upgrade the element array manually, starting from cntr2 and restart
          // the overlap check. At present it is only a warning to the user with
          // some relevant information.
          // Find the primitives and volumes for the overlapping elements
          // The element structure should also maintain information on the
          // volumes that an element belongs to.
          int prim1 = (EleArr + cntr1 - 1)->PrimitiveNb;
          int volele1 = VolRef1[prim1];
          int prim2 = (EleArr + cntr2 - 1)->PrimitiveNb;
          int volele2 = VolRef1[prim2];
 
          neBEMMessage("neBEMDiscretize - Overlapping collocation points!");
          printf("Element %d, primitive %d, volume %d overlaps with\n", cntr1,
                 prim1, volele1);
          printf("\telement %d, primitive %d, volume %d.\n", cntr2, prim2,
                 volele2);
          printf("\tposition 1: (%g , %g , %g) micron,\n", 1e6 * pt1.X,
                 1e6 * pt1.Y, 1e6 * pt1.Z);
          printf("\tposition 2: (%g , %g , %g) micron.\n", 1e6 * pt2.X,
                 1e6 * pt2.Y, 1e6 * pt2.Z);
          printf("Please redo the geometry.\n");
          return -1;
        }  // if dist <= MINDIST
      }    // for cntr2
    }      // for cntr1
  }        // check collocation point overlap
 
  NbElements = EleCntr;  // the final number of elements
  printf("Total final number of elements: %d\n", NbElements);
 
  // Store element related data in a file for a new mesh created, if opted for
  if (NewMesh && OptStoreElements) {
    if (OptFormattedFile) {
      int fstatus = WriteElements();
      if (fstatus) {
        neBEMMessage("neBEMDiscretize - problem writing Elements.\n");
        return -1;
      }
    }  // formatted file
 
    if (OptUnformattedFile) {
      neBEMMessage(
          "neBEMDiscretize - unformatted write not inplemented yet.\n");
      return -1;
    }  // unformatted file
  }    // store elements
 
  neBEMState = 4;
  stopClock = clock();
  neBEMTimeElapsed(startClock, stopClock);
  printf("to complete discretization\n");
 
  return (0);
}  // neBEMDiscretize ends

neBEMEnd()

INTFACEGLOBAL int neBEMEnd

(

void

)

Definition at line 2684 of file neBEMInterface.c.

                   {
  fprintf(fIsles,
          "IslesCntr: %d, ExactCntr: %d, FailureCntr: %d, ApproxCntr: %d\n",
          IslesCntr, ExactCntr, FailureCntr, ApproxCntr);
  fclose(fIsles);
  fIsles = NULL;
  printf("neBEM ends ... bye!\n");
 
  return 0;
}  // neBEMEnd ends

neBEMGetBoundingPlanes()

INTFACEGLOBAL int neBEMGetBoundingPlanes	(	int *	ixmin,
		double *	cxmin,
		double *	vxmin,
		int *	ixmax,
		double *	cxmax,
		double *	vxmax,
		int *	iymin,
		double *	cymin,
		double *	vymin,
		int *	iymax,
		double *	cymax,
		double *	vymax,
		int *	izmin,
		double *	czmin,
		double *	vzmin,
		int *	izmax,
		double *	czmax,
		double *	vzmax )

neBEMGetInputsFromFiles()

INTFACEGLOBAL int neBEMGetInputsFromFiles

(

void

)

neBEMGetMirror()

INTFACEGLOBAL int neBEMGetMirror	(	int	prim,
		int *	ix,
		int *	jx,
		double *	sx,
		int *	iy,
		int *	jy,
		double *	sy,
		int *	iz,
		int *	jz,
		double *	sz )

neBEMGetNbOfLines()

INTFACEGLOBAL int neBEMGetNbOfLines

(

const char *

fname

)

neBEMGetNbPrimitives()

INTFACEGLOBAL int neBEMGetNbPrimitives

(

void

)

neBEMGetPeriodicities()

INTFACEGLOBAL int neBEMGetPeriodicities	(	int	prim,
		int *	ix,
		int *	jx,
		double *	sx,
		int *	iy,
		int *	jy,
		double *	sy,
		int *	iz,
		int *	jz,
		double *	sz )

neBEMGetPrimitive()

INTFACEGLOBAL int neBEMGetPrimitive	(	int	prim,
		int *	nvertex,
		double	xvert[],
		double	yvert[],
		double	zvert[],
		double *	xnorm,
		double *	ynorm,
		double *	znorm,
		int *	volref1,
		int *	volref2 )

neBEMInitialize()

INTFACEGLOBAL int neBEMInitialize

(

void

)

Definition at line 39 of file neBEMInterface.c.

                          {
  // Version information
  strcpy(neBEMVersion, "1.9.09");
  strcpy(ISLESVersion, "1.4.9");
  printf("Using neBEM version %s and ISLES version %s\n", neBEMVersion,
         ISLESVersion);
 
  // The following is not an absolute necessity, but can ease the burden of the
  // user by setting default values of some of the important variable. Please
  // take a look at the src/Interface/ExampleDev2neBEM.c to find out what values
  // are being set, and what they mean.
  int fstatus = neBEMSetDefaults();
  if (fstatus != 0) {
    neBEMMessage("neBEMInitialize - neBEMSetDefaults");
    return -1;
  }
 
  // Change some of the global variables according to requirements.
  // Note that the solution flags and counters are set before the neBEMSolve
  // is invoked.
  LengthScale = 1.0;
  DebugLevel = 0;
 
  // The following function allows input files to be read and the geometry
  // set according to these files, as was customary in the previous versions
  if (OptDeviceFile) {
    printf("Reading geometry details from %s\n", DeviceInputFile);
    fstatus = neBEMGetInputsFromFiles();
    if (fstatus != 0) {
      neBEMMessage("neBEMInitialize - neBEMGetInputFromFiles");
      return -1;
    }
  }
 
  // creation of the mother output directory should be necessary only once
  if (neBEMState == 0) {
    fstatus = CreateDirStr();
    if (fstatus != 0) {
      neBEMMessage("neBEMInitialize - CreateDirStr");
      return -1;
    }
  }
 
  // Create Isles log file for keeping track of approximations (numerical
  // quadrature) in case evaluation of algebraic expressions fails.
  char IslesFile[256];
  strcpy(IslesFile, PPOutDir);
  strcat(IslesFile, "/Isles.log");
  fIsles = fopen(IslesFile, "w");
  if (fIsles == NULL) {
    neBEMMessage("neBEMInitialize - IslesFile");
    return -1;
  }
 
  // following integers keep track of the success / failure of the exact
  // expressions in estimating realistic physical properties.
  IslesCntr = ExactCntr = FailureCntr = ApproxCntr = 0;
 
  // Set up parameters related to neBEM computations
  // The file will be removed soon
#ifdef _OPENMP
  int RqstdThreads = 1;
  if (NbThreads > 0) RqstdThreads = NbThreads;
#endif
  /*
  FILE *processFile = fopen("neBEMInp/neBEMProcess.inp", "r");
  if (processFile == NULL) {
    printf("neBEMProcess.inp absent ... assuming defaults ...\n");
    PrimAfter = 0;
    if (NbThreads > 0) RqstdThreads = NbThreads;
  } else {
    fscanf(processFile, "PrimAfter: %d\n", &PrimAfter);
    fscanf(processFile, "RqstdThreads: %d\n", &RqstdThreads);
    fclose(processFile);
  }
  */
  // PrimAfter is assigned a value through "SetPrimAfter" member function of
  // ComponentNeBem3d class.
  // Default value of PrimAfter is a negative integer.
 
#ifdef _OPENMP
  int MaxProcessors = omp_get_num_procs();
  if (RqstdThreads > 1) {
    if (RqstdThreads < MaxProcessors) {
      // one processor left alone
      omp_set_num_threads(RqstdThreads);
    } else {
      printf("RqstdThreads: %d\n", RqstdThreads);
      RqstdThreads = MaxProcessors - 1;
      omp_set_num_threads(RqstdThreads);
      printf("Adjusted RqstdThreads: %d\n", RqstdThreads);
    }
  } else {
    // Work with one thread
    RqstdThreads = 1;  // cannot be zero or negative!
    omp_set_num_threads(RqstdThreads);
    printf("RqstdThreads: %d => No Multi-threading ...\n", RqstdThreads);
  }
 
  // OpenMP related information
  printf("PrimAfter: %d\n", PrimAfter);
  printf("RqstdThreads: %d, MaxProcessors: %d\n", RqstdThreads, MaxProcessors);
  printf("Maximum number of threads to be used for parallelization: %d\n",
         omp_get_max_threads());
  printf("Number of threads used for neBEMInitialize: %d\n",
         omp_get_num_threads());
#endif
 
  // Set up parameters related to voxelized data export for Garfield++
  // To be removed soon
  FILE *voxelInpFile = fopen("neBEMInp/neBEMVoxel.inp", "r");
  if (voxelInpFile == NULL) {
    printf("neBEMVoxel.inp absent ... assuming OptVoxel = 0 ...\n");
    OptVoxel = 0;
    OptStaggerVoxel = 0;
  } else {
    fscanf(voxelInpFile, "OptVoxel: %d\n", &OptVoxel);
    fscanf(voxelInpFile, "OptStaggerVoxel: %d\n", &OptStaggerVoxel);
    fscanf(voxelInpFile, "Xmin: %le\n", &Voxel.Xmin);
    fscanf(voxelInpFile, "Xmax: %le\n", &Voxel.Xmax);
    fscanf(voxelInpFile, "Ymin: %le\n", &Voxel.Ymin);
    fscanf(voxelInpFile, "Ymax: %le\n", &Voxel.Ymax);
    fscanf(voxelInpFile, "Zmin: %le\n", &Voxel.Zmin);
    fscanf(voxelInpFile, "Zmax: %le\n", &Voxel.Zmax);
    fscanf(voxelInpFile, "XStagger: %le\n", &Voxel.XStagger);
    fscanf(voxelInpFile, "YStagger: %le\n", &Voxel.YStagger);
    fscanf(voxelInpFile, "ZStagger: %le\n", &Voxel.ZStagger);
    fscanf(voxelInpFile, "NbOfXCells: %d\n", &Voxel.NbXCells);
    fscanf(voxelInpFile, "NbOfYCells: %d\n", &Voxel.NbYCells);
    fscanf(voxelInpFile, "NbOfZCells: %d\n", &Voxel.NbZCells);
    fclose(voxelInpFile);
  }  // inputs for Voxel
 
  // Set up parameters related to 3dMap data export for Garfield++
  // To be removed soon
  FILE *mapInpFile = fopen("neBEMInp/neBEMMap.inp", "r");
  if (mapInpFile == NULL) {
    printf("neBEMMap.inp absent ... assuming OptMap = 0 ...\n");
    OptMap = 0;
    OptStaggerMap = 0;
  } else {
    // While reading the input, OptMap and OptStaggerMap have to be read
    // first since that will decide whether there is a map and its version.
    fscanf(mapInpFile, "OptMap: %d\n", &OptMap);
    fscanf(mapInpFile, "OptStaggerMap: %d\n", &OptStaggerMap);
    fscanf(mapInpFile, "MapVersion: %9s\n", MapVersion);
    fscanf(mapInpFile, "Xmin: %le\n", &Map.Xmin);
    fscanf(mapInpFile, "Xmax: %le\n", &Map.Xmax);
    fscanf(mapInpFile, "Ymin: %le\n", &Map.Ymin);
    fscanf(mapInpFile, "Ymax: %le\n", &Map.Ymax);
    fscanf(mapInpFile, "Zmin: %le\n", &Map.Zmin);
    fscanf(mapInpFile, "Zmax: %le\n", &Map.Zmax);
    fscanf(mapInpFile, "XStagger: %le\n", &Map.XStagger);
    fscanf(mapInpFile, "YStagger: %le\n", &Map.YStagger);
    fscanf(mapInpFile, "ZStagger: %le\n", &Map.ZStagger);
    fscanf(mapInpFile, "NbOfXCells: %d\n", &Map.NbXCells);
    fscanf(mapInpFile, "NbOfYCells: %d\n", &Map.NbYCells);
    fscanf(mapInpFile, "NbOfZCells: %d\n", &Map.NbZCells);
    fclose(mapInpFile);
  }  // inputs for 3dMap
 
  // Set up parameters related to fast volume
  if (OptFastVol) {
    FILE *fastInpFile = fopen("neBEMInp/neBEMFastVol.inp", "r");
    if (fastInpFile == NULL) {
      printf("neBEMFastVol.inp absent ... assuming OptFastVol = 0 ...\n");
      OptFastVol = 0;
      OptStaggerFastVol = 0;
      OptCreateFastPF = 0;
      OptReadFastPF = 0;
      FastVol.NbBlocks = 0;
      FastVol.NbOmitVols = 0;
      FastVol.NbIgnoreVols = 0;
    } else {
      fscanf(fastInpFile, "OptFastVol: %d\n", &OptFastVol);
      fscanf(fastInpFile, "OptStaggerFastVol: %d\n", &OptStaggerFastVol);
      fscanf(fastInpFile, "OptCreateFastPF: %d\n", &OptCreateFastPF);
      fscanf(fastInpFile, "OptReadFastPF: %d\n", &OptReadFastPF);
      fscanf(fastInpFile, "NbPtSkip: %d\n", &NbPtSkip);
      fscanf(fastInpFile, "NbStgPtSkip: %d\n", &NbStgPtSkip);
      fscanf(fastInpFile, "LX: %le\n", &FastVol.LX);
      fscanf(fastInpFile, "LY: %le\n", &FastVol.LY);
      fscanf(fastInpFile, "LZ: %le\n", &FastVol.LZ);
      fscanf(fastInpFile, "CornerX: %le\n", &FastVol.CrnrX);
      fscanf(fastInpFile, "CornerY: %le\n", &FastVol.CrnrY);
      fscanf(fastInpFile, "CornerZ: %le\n", &FastVol.CrnrZ);
      fscanf(fastInpFile, "YStagger: %le\n", &FastVol.YStagger);
      if (!OptStaggerFastVol)
        FastVol.YStagger = 0.0;  // ignore any non-zero value
      fscanf(fastInpFile, "NbOfBlocks: %d\n", &FastVol.NbBlocks);
      BlkNbXCells = ivector(1, FastVol.NbBlocks);
      BlkNbYCells = ivector(1, FastVol.NbBlocks);
      BlkNbZCells = ivector(1, FastVol.NbBlocks);
      BlkLZ = dvector(1, FastVol.NbBlocks);
      BlkCrnrZ = dvector(1, FastVol.NbBlocks);
      for (int block = 1; block <= FastVol.NbBlocks; ++block) {
        fscanf(fastInpFile, "NbOfXCells: %d\n", &BlkNbXCells[block]);
        fscanf(fastInpFile, "NbOfYCells: %d\n", &BlkNbYCells[block]);
        fscanf(fastInpFile, "NbOfZCells: %d\n", &BlkNbZCells[block]);
        fscanf(fastInpFile, "LZ: %le\n", &BlkLZ[block]);
        fscanf(fastInpFile, "CornerZ: %le\n", &BlkCrnrZ[block]);
      }  // inputs for blocks
      fscanf(fastInpFile, "NbOfOmitVols: %d\n", &FastVol.NbOmitVols);
      if (FastVol.NbOmitVols) {
        OmitVolLX = dvector(1, FastVol.NbOmitVols);
        OmitVolLY = dvector(1, FastVol.NbOmitVols);
        OmitVolLZ = dvector(1, FastVol.NbOmitVols);
        OmitVolCrnrX = dvector(1, FastVol.NbOmitVols);
        OmitVolCrnrY = dvector(1, FastVol.NbOmitVols);
        OmitVolCrnrZ = dvector(1, FastVol.NbOmitVols);
        for (int omit = 1; omit <= FastVol.NbOmitVols; ++omit) {
          fscanf(fastInpFile, "OmitVolLX: %le\n", &OmitVolLX[omit]);
          fscanf(fastInpFile, "OmitVolLY: %le\n", &OmitVolLY[omit]);
          fscanf(fastInpFile, "OmitVolLZ: %le\n", &OmitVolLZ[omit]);
          fscanf(fastInpFile, "OmitVolCornerX: %le\n", &OmitVolCrnrX[omit]);
          fscanf(fastInpFile, "OmitVolCornerY: %le\n", &OmitVolCrnrY[omit]);
          fscanf(fastInpFile, "OmitVolCornerZ: %le\n", &OmitVolCrnrZ[omit]);
        }  // inputs for OmitVols
      }    // inputs for OmitVols
      fscanf(fastInpFile, "NbOfIgnoreVols: %d\n", &FastVol.NbIgnoreVols);
      if (FastVol.NbIgnoreVols) {
        IgnoreVolLX = dvector(1, FastVol.NbIgnoreVols);
        IgnoreVolLY = dvector(1, FastVol.NbIgnoreVols);
        IgnoreVolLZ = dvector(1, FastVol.NbIgnoreVols);
        IgnoreVolCrnrX = dvector(1, FastVol.NbIgnoreVols);
        IgnoreVolCrnrY = dvector(1, FastVol.NbIgnoreVols);
        IgnoreVolCrnrZ = dvector(1, FastVol.NbIgnoreVols);
        for (int ignore = 1; ignore <= FastVol.NbIgnoreVols; ++ignore) {
          fscanf(fastInpFile, "IgnoreVolLX: %le\n", &IgnoreVolLX[ignore]);
          fscanf(fastInpFile, "IgnoreVolLY: %le\n", &IgnoreVolLY[ignore]);
          fscanf(fastInpFile, "IgnoreVolLZ: %le\n", &IgnoreVolLZ[ignore]);
          fscanf(fastInpFile, "IgnoreVolCornerX: %le\n",
                 &IgnoreVolCrnrX[ignore]);
          fscanf(fastInpFile, "IgnoreVolCornerY: %le\n",
                 &IgnoreVolCrnrY[ignore]);
          fscanf(fastInpFile, "IgnoreVolCornerZ: %le\n",
                 &IgnoreVolCrnrZ[ignore]);
        }  // inputs for IgnoreVols
      }    // inputs for IgnoreVols
      for (int ignore = 1; ignore <= FastVol.NbIgnoreVols; ++ignore) {
        printf("IgnoreVolLX: %le\n", IgnoreVolLX[ignore]);
        printf("IgnoreVolLY: %le\n", IgnoreVolLY[ignore]);
        printf("IgnoreVolLZ: %le\n", IgnoreVolLZ[ignore]);
        printf("IgnoreVolCornerX: %le\n", IgnoreVolCrnrX[ignore]);
        printf("IgnoreVolCornerY: %le\n", IgnoreVolCrnrY[ignore]);
        printf("IgnoreVolCornerZ: %le\n", IgnoreVolCrnrZ[ignore]);
      }  // inputs for IgnoreVols
      fclose(fastInpFile);
    }  // else fastInpFile
  }    // if OptFastVol
 
  printf("neBEM initialized ...\n");
  fflush(stdout);
  sleep(3);  // wait for three seconds so that the user gets time to react
 
  neBEMState = 1;  // state 1 implied initialization of neBEM completed
 
  // announce success - later, add the name of the calling code
  return (0);
}  // neBEMInitialize ends

neBEMMessage()

INTFACEGLOBAL int neBEMMessage

(

const char *

message

)

Definition at line 2855 of file neBEMInterface.c.

                                      {
  fprintf(stdout, "neBEMMessage: %s\n", message);
 
  return 0;
}  // neBEMMessage ends

Referenced by AnalyzePrimitive(), ComputeSolution(), CreateDirStr(), CreateFastVolElePF(), CreateWtFldFastVolPF(), DiscretizeRectangle(), DiscretizeTriangle(), DiscretizeWire(), FastKnChPFAtPoint(), FastPFAtPoint(), InitChargingUp(), InitialConditions(), InitKnownCharges(), InvertMatrix(), LHMatrix(), MapFPR(), neBEMBoundaryInitialConditions(), neBEMDiscretize(), neBEMInitialize(), neBEMPF(), neBEMPrepareWeightingField(), neBEMReadGeometry(), neBEMSolve(), neBEMWeightingField(), ReadElements(), neBEM::ReadInitFile(), ReadInvertedMatrix(), ReadPrimitives(), ReadSolution(), Solve(), SurfaceElements(), VoxelFPR(), WireElements(), WriteElements(), WritePrimitives(), and WtFldFastPFAtPoint().

neBEMPF()

INTFACEGLOBAL int neBEMPF	(	Point3D *	point,
		double *	potential,
		Vector3D *	field )

Definition at line 2096 of file neBEMInterface.c.

                                                                {
  if (neBEMState < 9) {
    printf("neBEMPF cannot be called before reaching state 9.\n");
    return (-1);
  }
 
  // printf("neBEMPF called %8d times", ++neBEMPFCallCntr);
 
  double Pot;
  int fstatus;
  if (OptFastVol)  // Note: this is not the Create or Read option
  {
    fstatus = FastPFAtPoint(point, &Pot, field);
    if (fstatus != 0) {
      neBEMMessage("neBEMPF - FastPFAtPoint");
      return -1;
    }
  } else {
    fstatus = PFAtPoint(point, &Pot, field);
    if (fstatus != 0) {
      neBEMMessage("neBEMPF - PFAtPoint");
      return -1;
    }
  }
 
  *potential = Pot;
 
  // printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
 
  return (0);
}  // neBEMPF ends

neBEMPrepareWeightingField()

INTFACEGLOBAL int neBEMPrepareWeightingField	(	int	NbPrimsWtField,
		int	PrimListWtField[] )

Definition at line 2141 of file neBEMInterface.c.

                                                          {
  static int IdWtField = 0;
 
  int dbgFn = 0;
  int fstatus = 0;
 
  if (neBEMState < 7) {
    printf(
        "neBEMPrepareWeightingField: Weighting computations only meaningful "
        "beyond neBEMState 7 ...\n");
    return -1;
  }
 
  // Find first free slot
  const int MaxWtField =
      MAXWtFld;  // used also while deallocating these memories
  if (WtFieldChDen == NULL)
    WtFieldChDen = (double **)malloc(MaxWtField * sizeof(double *));
  if (AvWtChDen == NULL)
    AvWtChDen = (double **)malloc(MaxWtField * sizeof(double *));
 
  ++IdWtField;
  if (IdWtField >= MaxWtField) {
    printf(
        "neBEMPrepareWeightingField: reached MaxWtField (%d) weighting "
        "fields.\n",
        MAXWtFld);
    return -1;
  } else {
    printf("\nPreparing weighting field for %d-th set.\n", IdWtField);
  }  // else within MaxWtField
 
  // Allocate a new column to store this solution set
  WtFieldChDen[IdWtField] = (double *)malloc((NbElements + 2) * sizeof(double));
  AvWtChDen[IdWtField] = (double *)malloc((NbPrimitives + 2) * sizeof(double));
 
  fstatus = WeightingFieldSolution(nprim, primlist, WtFieldChDen[IdWtField]);
  if (fstatus) {
    neBEMMessage("neBEMPrepareWeightingField - WeightingFieldSolution");
    return -1;
  } else {
    printf("Computed weighting field solution\n");
  }
 
  // estimate primitive related avrg wt field charge densities
  // OMPCheck - may be parallelized
  for (int prim = 1; prim <= NbPrimitives; ++prim) {
    double area = 0.0;  // need area of the primitive as well!
    AvWtChDen[IdWtField][prim] = 0.0;
 
    for (int ele = ElementBgn[prim]; ele <= ElementEnd[prim]; ++ele) {
      area += (EleArr + ele - 1)->G.dA;
      AvWtChDen[IdWtField][prim] +=
          WtFieldChDen[IdWtField][ele] * (EleArr + ele - 1)->G.dA;
    }
 
    AvWtChDen[IdWtField][prim] /= area;
  }
  printf("Computed primitive-averaged weighting field solutions\n");
 
  // stringify the integer
  char strIdWtField[5];
  snprintf(strIdWtField, 5, "%d", IdWtField);
  // printf("strIdWtField: %s\n", strIdWtField);
 
  // Set up parameters related to fixed specification of weighting field
  // printf("OptFixedWtField: %d\n", OptFixedWtField[IdWtField]);fflush(stdout);
  if (OptFixedWtField[IdWtField]) {
    char fileWtField[256];
    strcpy(fileWtField, "neBEMInp/neBEMFixedWtField_");
    strcat(fileWtField, strIdWtField);
    strcat(fileWtField, ".inp");
    FILE *fixedWtInpFile = fopen(fileWtField, "r");
    if (fixedWtInpFile == NULL) {
      printf(
          "neBEMFixedWtField.inp absent ... assuming OptFixedWtField = 0 "
          "...\n");
      OptFixedWtField[IdWtField] = 0;
      FixedWtPotential[IdWtField] = 0.0;
      FixedWtFieldX[IdWtField] = 0.0;
      FixedWtFieldY[IdWtField] = 0.0;
      FixedWtFieldZ[IdWtField] = 0.0;
    } else {
      fscanf(fixedWtInpFile, "OptFixedWtField: %d\n",
             &OptFixedWtField[IdWtField]);
      fscanf(fixedWtInpFile, "FixedWtPotential: %lg\n",
             &FixedWtPotential[IdWtField]);
      fscanf(fixedWtInpFile, "FixedWtFieldX: %lg\n", &FixedWtFieldX[IdWtField]);
      fscanf(fixedWtInpFile, "FixedWtFieldY: %lg\n", &FixedWtFieldY[IdWtField]);
      fscanf(fixedWtInpFile, "FixedWtFieldZ: %lg\n", &FixedWtFieldZ[IdWtField]);
      fclose(fixedWtInpFile);
    }  // else fixedWtInpFile
  }    // if OptFixedWtField
 
  // Weighting field fast volume related computations
  // Set up parameters related to weighting field fast volumes.
  // There can be multiple weighting fields depending on readout electrodes.
  // As a result, there can be multiple versions of this input file, each
  // reflecting the requirement of one set of electrodes.
  // printf("OptWtFldFastVol: %d\n", OptWtFldFastVol[IdWtField]);fflush(stdout);
  if (OptWtFldFastVol[IdWtField]) {
    char fileWtField[256];
    strcpy(fileWtField, "neBEMInp/neBEMWtFldFastVol_");
    strcat(fileWtField, strIdWtField);
    strcat(fileWtField, ".inp");
    FILE *fastWtFldInpFile = fopen(fileWtField, "r");
    // FILE *fastWtFldInpFile = fopen("neBEMInp/neBEMWtFldFastVol.inp", "r");
    if (fastWtFldInpFile == NULL) {
      printf(
          "neBEMWtFldFastVol.inp absent ... assuming OptWtFldFastVol = 0 "
          "...\n");
      OptWtFldFastVol[IdWtField] = 0;
      OptStaggerWtFldFastVol[IdWtField] = 0;
      OptCreateWtFldFastPF[IdWtField] = 0;
      OptReadWtFldFastPF[IdWtField] = 0;
      WtFldFastVol[IdWtField].NbBlocks = 0;
      WtFldFastVol[IdWtField].NbOmitVols = 0;
      WtFldFastVol[IdWtField].NbIgnoreVols = 0;
    } else {
      fscanf(fastWtFldInpFile, "OptFastVol: %d\n", &OptWtFldFastVol[IdWtField]);
      fscanf(fastWtFldInpFile, "OptStaggerFastVol: %d\n",
             &OptStaggerWtFldFastVol[IdWtField]);
      fscanf(fastWtFldInpFile, "OptCreateFastPF: %d\n",
             &OptCreateWtFldFastPF[IdWtField]);
      fscanf(fastWtFldInpFile, "OptReadFastPF: %d\n",
             &OptReadWtFldFastPF[IdWtField]);
      fscanf(fastWtFldInpFile, "NbPtSkip: %d\n", &WtFldNbPtSkip[IdWtField]);
      fscanf(fastWtFldInpFile, "NbStgPtSkip: %d\n",
             &StgWtFldNbPtSkip[IdWtField]);
      fscanf(fastWtFldInpFile, "LX: %le\n", &WtFldFastVol[IdWtField].LX);
      fscanf(fastWtFldInpFile, "LY: %le\n", &WtFldFastVol[IdWtField].LY);
      fscanf(fastWtFldInpFile, "LZ: %le\n", &WtFldFastVol[IdWtField].LZ);
      fscanf(fastWtFldInpFile, "CornerX: %le\n",
             &WtFldFastVol[IdWtField].CrnrX);
      fscanf(fastWtFldInpFile, "CornerY: %le\n",
             &WtFldFastVol[IdWtField].CrnrY);
      fscanf(fastWtFldInpFile, "CornerZ: %le\n",
             &WtFldFastVol[IdWtField].CrnrZ);
      fscanf(fastWtFldInpFile, "YStagger: %le\n",
             &WtFldFastVol[IdWtField].YStagger);
      if (!OptStaggerWtFldFastVol[IdWtField])
        WtFldFastVol[IdWtField].YStagger = 0.0;  // ignore any non-zero value
      fscanf(fastWtFldInpFile, "NbOfBlocks: %d\n",
             &WtFldFastVol[IdWtField].NbBlocks);
      WtFldBlkNbXCells[IdWtField] =
          ivector(1, WtFldFastVol[IdWtField].NbBlocks);
      WtFldBlkNbYCells[IdWtField] =
          ivector(1, WtFldFastVol[IdWtField].NbBlocks);
      WtFldBlkNbZCells[IdWtField] =
          ivector(1, WtFldFastVol[IdWtField].NbBlocks);
      WtFldBlkLZ[IdWtField] = dvector(1, WtFldFastVol[IdWtField].NbBlocks);
      WtFldBlkCrnrZ[IdWtField] = dvector(1, WtFldFastVol[IdWtField].NbBlocks);
      for (int block = 1; block <= WtFldFastVol[IdWtField].NbBlocks; ++block) {
        fscanf(fastWtFldInpFile, "NbOfXCells: %d\n",
               &WtFldBlkNbXCells[IdWtField][block]);
        fscanf(fastWtFldInpFile, "NbOfYCells: %d\n",
               &WtFldBlkNbYCells[IdWtField][block]);
        fscanf(fastWtFldInpFile, "NbOfZCells: %d\n",
               &WtFldBlkNbZCells[IdWtField][block]);
        fscanf(fastWtFldInpFile, "LZ: %le\n", &WtFldBlkLZ[IdWtField][block]);
        fscanf(fastWtFldInpFile, "CornerZ: %le\n",
               &WtFldBlkCrnrZ[IdWtField][block]);
      }  // inputs for blocks
      fscanf(fastWtFldInpFile, "NbOfOmitVols: %d\n",
             &WtFldFastVol[IdWtField].NbOmitVols);
      if (WtFldFastVol[IdWtField].NbOmitVols) {
        WtFldOmitVolLX[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbOmitVols);
        WtFldOmitVolLY[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbOmitVols);
        WtFldOmitVolLZ[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbOmitVols);
        WtFldOmitVolCrnrX[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbOmitVols);
        WtFldOmitVolCrnrY[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbOmitVols);
        WtFldOmitVolCrnrZ[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbOmitVols);
        for (int omit = 1; omit <= WtFldFastVol[IdWtField].NbOmitVols; ++omit) {
          fscanf(fastWtFldInpFile, "OmitVolLX: %le\n",
                 &WtFldOmitVolLX[IdWtField][omit]);
          fscanf(fastWtFldInpFile, "OmitVolLY: %le\n",
                 &WtFldOmitVolLY[IdWtField][omit]);
          fscanf(fastWtFldInpFile, "OmitVolLZ: %le\n",
                 &WtFldOmitVolLZ[IdWtField][omit]);
          fscanf(fastWtFldInpFile, "OmitVolCornerX: %le\n",
                 &WtFldOmitVolCrnrX[IdWtField][omit]);
          fscanf(fastWtFldInpFile, "OmitVolCornerY: %le\n",
                 &WtFldOmitVolCrnrY[IdWtField][omit]);
          fscanf(fastWtFldInpFile, "OmitVolCornerZ: %le\n",
                 &WtFldOmitVolCrnrZ[IdWtField][omit]);
        }  // for loop inputs for OmitVols
      }    // inputs for OmitVols
      fscanf(fastWtFldInpFile, "NbOfIgnoreVols: %d\n",
             &WtFldFastVol[IdWtField].NbIgnoreVols);
      if (WtFldFastVol[IdWtField].NbIgnoreVols) {
        WtFldIgnoreVolLX[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbIgnoreVols);
        WtFldIgnoreVolLY[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbIgnoreVols);
        WtFldIgnoreVolLZ[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbIgnoreVols);
        WtFldIgnoreVolCrnrX[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbIgnoreVols);
        WtFldIgnoreVolCrnrY[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbIgnoreVols);
        WtFldIgnoreVolCrnrZ[IdWtField] =
            dvector(1, WtFldFastVol[IdWtField].NbIgnoreVols);
        for (int ignore = 1; ignore <= WtFldFastVol[IdWtField].NbIgnoreVols;
             ++ignore) {
          fscanf(fastWtFldInpFile, "IgnoreVolLX: %le\n",
                 &WtFldIgnoreVolLX[IdWtField][ignore]);
          fscanf(fastWtFldInpFile, "IgnoreVolLY: %le\n",
                 &WtFldIgnoreVolLY[IdWtField][ignore]);
          fscanf(fastWtFldInpFile, "IgnoreVolLZ: %le\n",
                 &WtFldIgnoreVolLZ[IdWtField][ignore]);
          fscanf(fastWtFldInpFile, "IgnoreVolCornerX: %le\n",
                 &WtFldIgnoreVolCrnrX[IdWtField][ignore]);
          fscanf(fastWtFldInpFile, "IgnoreVolCornerY: %le\n",
                 &WtFldIgnoreVolCrnrY[IdWtField][ignore]);
          fscanf(fastWtFldInpFile, "IgnoreVolCornerZ: %le\n",
                 &WtFldIgnoreVolCrnrZ[IdWtField][ignore]);
        }  // for loop inputs for IgnoreVols
      }    // inputs for IgnoreVols
      if (dbgFn) {
        for (int ignore = 1; ignore <= WtFldFastVol[IdWtField].NbIgnoreVols;
             ++ignore) {
          printf("WtFldIgnoreVolLX: %le\n",
                 WtFldIgnoreVolLX[IdWtField][ignore]);
          printf("WtFldIgnoreVolLY: %le\n",
                 WtFldIgnoreVolLY[IdWtField][ignore]);
          printf("WtFldIgnoreVolLZ: %le\n",
                 WtFldIgnoreVolLZ[IdWtField][ignore]);
          printf("WtFldIgnoreVolCornerX: %le\n",
                 WtFldIgnoreVolCrnrX[IdWtField][ignore]);
          printf("WtFldIgnoreVolCornerY: %le\n",
                 WtFldIgnoreVolCrnrY[IdWtField][ignore]);
          printf("WtFldIgnoreVolCornerZ: %le\n",
                 WtFldIgnoreVolCrnrZ[IdWtField][ignore]);
        }  // inputs for IgnoreVols
      }
      fclose(fastWtFldInpFile);
    }  // else fastWtFldInpFile
 
    int MaxXCells = WtFldBlkNbXCells[IdWtField][1];
    int MaxYCells = WtFldBlkNbYCells[IdWtField][1];
    int MaxZCells = WtFldBlkNbZCells[IdWtField][1];
    clock_t startFastClock = clock();
 
    // find maximum number of Xcells etc in all the blocks
    // simplifies memory allocation using nrutils but hogs memory!
    for (int block = 1; block <= WtFldFastVol[IdWtField].NbBlocks; ++block) {
      if (block == 1) {
        MaxXCells = WtFldBlkNbXCells[IdWtField][1];
        MaxYCells = WtFldBlkNbYCells[IdWtField][1];
        MaxZCells = WtFldBlkNbZCells[IdWtField][1];
      } else {
        if (MaxXCells < WtFldBlkNbXCells[IdWtField][block])
          MaxXCells = WtFldBlkNbXCells[IdWtField][block];
        if (MaxYCells < WtFldBlkNbYCells[IdWtField][block])
          MaxYCells = WtFldBlkNbYCells[IdWtField][block];
        if (MaxZCells < WtFldBlkNbZCells[IdWtField][block])
          MaxZCells = WtFldBlkNbZCells[IdWtField][block];
      }
    }  // loop block for finding maxm cells among all the blocks
 
    if (dbgFn) {
      printf("OptWtFldFastVol: %d\n", OptWtFldFastVol[IdWtField]);
      printf("OptStaggerWtFldFastVol: %d\n", OptStaggerWtFldFastVol[IdWtField]);
      printf("OptCreateWtFldFastPF: %d\n", OptCreateWtFldFastPF[IdWtField]);
      printf("OptReadWtFldFastPF: %d\n", OptReadWtFldFastPF[IdWtField]);
      printf("WtFldNbPtSkip: %d\n", WtFldNbPtSkip[IdWtField]);
      printf("StgWtFldNbPtSkip: %d\n", StgWtFldNbPtSkip[IdWtField]);
      printf("LX: %le\n", WtFldFastVol[IdWtField].LX);
      printf("LY: %le\n", WtFldFastVol[IdWtField].LY);
      printf("LZ: %le\n", WtFldFastVol[IdWtField].LZ);
      printf("CornerX: %le\n", WtFldFastVol[IdWtField].CrnrX);
      printf("CornerY: %le\n", WtFldFastVol[IdWtField].CrnrY);
      printf("CornerZ: %le\n", WtFldFastVol[IdWtField].CrnrZ);
      printf("YStagger: %le\n", WtFldFastVol[IdWtField].YStagger);
      printf("NbOfBlocks: %d\n", WtFldFastVol[IdWtField].NbBlocks);
      for (int block = 1; block <= WtFldFastVol[IdWtField].NbBlocks; ++block) {
        printf("NbOfXCells[%d]: %d\n", block,
               WtFldBlkNbXCells[IdWtField][block]);
        printf("NbOfYCells[%d]: %d\n", block,
               WtFldBlkNbYCells[IdWtField][block]);
        printf("NbOfZCells[%d]: %d\n", block,
               WtFldBlkNbZCells[IdWtField][block]);
        printf("LZ[%d]: %le\n", block, WtFldBlkLZ[IdWtField][block]);
        printf("CornerZ[%d]: %le\n", block, WtFldBlkCrnrZ[IdWtField][block]);
      }
      printf("NbOfOmitVols: %d\n", WtFldFastVol[IdWtField].NbOmitVols);
      if (WtFldFastVol[IdWtField].NbOmitVols) {
        for (int omit = 1; omit <= WtFldFastVol[IdWtField].NbOmitVols; ++omit) {
          printf("OmitVolLX[%d]: %le\n", omit, WtFldOmitVolLX[IdWtField][omit]);
          printf("OmitVolLY[%d]: %le\n", omit, WtFldOmitVolLY[IdWtField][omit]);
          printf("OmitVolLZ[%d]: %le\n", omit, WtFldOmitVolLZ[IdWtField][omit]);
          printf("OmitCrnrX[%d]: %le\n", omit,
                 WtFldOmitVolCrnrX[IdWtField][omit]);
          printf("OmitCrnrY[%d]: %le\n", omit,
                 WtFldOmitVolCrnrY[IdWtField][omit]);
          printf("OmitCrnrZ[%d]: %le\n", omit,
                 WtFldOmitVolCrnrZ[IdWtField][omit]);
        }
      }
      printf("MaxXCells, MaxYCells, MaxZCells: %d, %d, %d\n", MaxXCells,
             MaxYCells, MaxZCells);
      printf("NbOfIgnoreVols: %d\n", WtFldFastVol[IdWtField].NbIgnoreVols);
    }  // dbgFn
 
    // Following memory allocations are necessary even for creating
    // the fast volumes.
    // In fact, since these tensors are already in place, the fast volumes
    // can be used to evaluate potential and fields immediately after they are
    // created.
    /* Memory wastage!!! Improve as soon as possible. */
    WtFldFastPot[IdWtField] =
        d4tensor(1, WtFldFastVol[IdWtField].NbBlocks, 1, MaxXCells + 1, 1,
                 MaxYCells + 1, 1, MaxZCells + 1);
    WtFldFastFX[IdWtField] =
        d4tensor(1, WtFldFastVol[IdWtField].NbBlocks, 1, MaxXCells + 1, 1,
                 MaxYCells + 1, 1, MaxZCells + 1);
    WtFldFastFY[IdWtField] =
        d4tensor(1, WtFldFastVol[IdWtField].NbBlocks, 1, MaxXCells + 1, 1,
                 MaxYCells + 1, 1, MaxZCells + 1);
    WtFldFastFZ[IdWtField] =
        d4tensor(1, WtFldFastVol[IdWtField].NbBlocks, 1, MaxXCells + 1, 1,
                 MaxYCells + 1, 1, MaxZCells + 1);
 
    if (OptStaggerWtFldFastVol[IdWtField]) {
      /* Memory wastage!!! Improve as soon as possible. */
      StgWtFldFastPot[IdWtField] =
          d4tensor(1, WtFldFastVol[IdWtField].NbBlocks, 1, MaxXCells + 1, 1,
                   MaxYCells + 1, 1, MaxZCells + 1);
      StgWtFldFastFX[IdWtField] =
          d4tensor(1, WtFldFastVol[IdWtField].NbBlocks, 1, MaxXCells + 1, 1,
                   MaxYCells + 1, 1, MaxZCells + 1);
      StgWtFldFastFY[IdWtField] =
          d4tensor(1, WtFldFastVol[IdWtField].NbBlocks, 1, MaxXCells + 1, 1,
                   MaxYCells + 1, 1, MaxZCells + 1);
      StgWtFldFastFZ[IdWtField] =
          d4tensor(1, WtFldFastVol[IdWtField].NbBlocks, 1, MaxXCells + 1, 1,
                   MaxYCells + 1, 1, MaxZCells + 1);
    }  // if OptStaggerWtFldFastVol
 
    if (OptCreateWtFldFastPF[IdWtField]) {
      fstatus = CreateWtFldFastVolPF(IdWtField);
 
      clock_t stopFastClock = clock();
      neBEMTimeElapsed(startFastClock, stopFastClock);
      printf("to compute WtFldFastVolPF\n");
      // neBEMMessage(
      // "neBEMSolve - Failure computing WtFldFastVolPF: not implemented");
      // return -1;
    }  // if OptCreateWtFldFastPF
 
    if (OptReadWtFldFastPF[IdWtField]) {  // reading from file
      int nbXCells, nbYCells, nbZCells;
      int tmpblk;
      double xpt, ypt, zpt;
 
      // stringify the integer
      char stringIdWtField[5];
      snprintf(stringIdWtField, 5, "%d", IdWtField);
      char FastVolPFFile[256];
      strcpy(FastVolPFFile, BCOutDir);
      strcat(FastVolPFFile, "/WtFldFastVolPF_");
      strcat(FastVolPFFile, stringIdWtField);
      strcat(FastVolPFFile, ".out");
      FILE *fFastVolPF = fopen(FastVolPFFile, "r");
      if (fFastVolPF == NULL) {
        neBEMMessage("in neBEMSolve - WtFldFastVolPFFile");
        return -1;
      }
 
      fscanf(fFastVolPF, "#block\tX\tY\tZ\tPot\tFX\tFY\tFZ\n");
 
      for (int block = 1; block <= WtFldFastVol[IdWtField].NbBlocks; ++block) {
        nbXCells = WtFldBlkNbXCells[IdWtField][block];
        nbYCells = WtFldBlkNbYCells[IdWtField][block];
        nbZCells = WtFldBlkNbZCells[IdWtField][block];
 
        for (int i = 1; i <= nbXCells + 1; ++i) {
          for (int j = 1; j <= nbYCells + 1; ++j) {
            for (int k = 1; k <= nbZCells + 1; ++k) {
              fscanf(fFastVolPF, "%d\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",
                     &tmpblk, &xpt, &ypt, &zpt,
                     &WtFldFastPot[IdWtField][block][i][j][k],
                     &WtFldFastFX[IdWtField][block][i][j][k],
                     &WtFldFastFY[IdWtField][block][i][j][k],
                     &WtFldFastFZ[IdWtField][block][i][j][k]);
            }  // loop k
          }    // loop j
        }      // loop i
      }        // loop block
      fclose(fFastVolPF);
 
      if (OptStaggerWtFldFastVol[IdWtField]) {
        // stringify the integer
        snprintf(stringIdWtField, 5, "%d", IdWtField);
        char StgFastVolPFFile[256];
        strcpy(StgFastVolPFFile, BCOutDir);
        // strcat(StgFastVolPFFile, "/WtFldStgFastVolPF.out");
        strcat(StgFastVolPFFile, "/StgWtFldFastVolPF_");
        strcat(StgFastVolPFFile, stringIdWtField);
        strcat(StgFastVolPFFile, ".out");
        FILE* fStgFastVolPF = fopen(StgFastVolPFFile, "r");
        if (fStgFastVolPF == NULL) {
          neBEMMessage("in neBEMSolve - StgWtFldFastVolPFFile");
          return -1;
        }
 
        fscanf(fStgFastVolPF, "#block\tX\tY\tZ\tPot\tFX\tFY\tFZ\n");
 
        for (int block = 1; block <= WtFldFastVol[IdWtField].NbBlocks;
             ++block) {
          nbXCells = WtFldBlkNbXCells[IdWtField][block];
          nbYCells = WtFldBlkNbYCells[IdWtField][block];
          nbZCells = WtFldBlkNbZCells[IdWtField][block];
 
          for (int i = 1; i <= nbXCells + 1; ++i) {
            for (int j = 1; j <= nbYCells + 1; ++j) {
              for (int k = 1; k <= nbZCells + 1; ++k) {
                fscanf(fStgFastVolPF, "%d\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",
                       &tmpblk, &xpt, &ypt, &zpt,
                       &StgWtFldFastPot[IdWtField][block][i][j][k],
                       &StgWtFldFastFX[IdWtField][block][i][j][k],
                       &StgWtFldFastFY[IdWtField][block][i][j][k],
                       &StgWtFldFastFZ[IdWtField][block][i][j][k]);
              }  // loop k
            }    // loop j
          }      // loop i
        }        // loop block
        fclose(fStgFastVolPF);
      }  // if OptStaggerWtFldFastVol
 
      clock_t stopFastClock = clock();
      neBEMTimeElapsed(startFastClock, stopFastClock);
      printf("to read WtFldFastVolPF\n");
    }  // if OptReadWtFldFastPF
  }    // if OptWtFldFastVol
 
  return IdWtField;
}  // neBEMPrepareWeightingField ends

neBEMReadGeometry()

INTFACEGLOBAL int neBEMReadGeometry

(

void

)

Definition at line 307 of file neBEMInterface.c.

                            {
  int dbgFn = 0;
  int fstatus;
 
  startClock = clock();
 
  // For a model that was defined before and for which data was stored in a file
  if ((!NewModel) && (!NewBC) && (OptStorePrimitives)) {
    fstatus = ReadPrimitives();
    if (fstatus) {
      neBEMMessage("neBEMReadGeometry - problem reading stored Primitives.\n");
      return -1;
    }
    neBEMState = 3;  // primitives read in after initialization and Nbs
    return 0;
  }
 
  printf("geometry inputs ...\n");
  if (neBEMState != 1) {
    printf("reading geometry possible only after initialization ...\n");
    return -1;
  }
 
  NbPrimitives = neBEMGetNbPrimitives();
  OrgnlNbPrimitives = NbPrimitives;
  if (NbPrimitives == 0) {
    // nothing to do - return control to calling routine
    neBEMMessage("neBEMReadGeometry - no primitive.\n");
    return (-1);  // for the time being
  }
 
  NbSurfs = 0;
  NbWires = 0;
  neBEMState = 2;
 
  // Allocate memory for storing the geometry primitives till the elements are
  // created.
  // Explicit storage of these variables related to primitives may not be
  // necessary if the elements are created immediately after a primitive is read
  // in.
  // MaxNbVertices = 4; // value specified through SetDefaults or init files
 
  // neBEM has been initialized, NbPrimitives set
  PrimType = ivector(1, NbPrimitives);
  NbVertices = ivector(1, NbPrimitives);
  OrgnlToEffPrim = imatrix(1, NbPrimitives, 0, 2);  // 0 init, 1 intrfc, 2 rmv
  XVertex = dmatrix(1, NbPrimitives, 0, MaxNbVertices - 1);
  YVertex = dmatrix(1, NbPrimitives, 0, MaxNbVertices - 1);
  ZVertex = dmatrix(1, NbPrimitives, 0, MaxNbVertices - 1);
  XNorm = dvector(1, NbPrimitives);
  YNorm = dvector(1, NbPrimitives);
  ZNorm = dvector(1, NbPrimitives);
  PrimLX = dvector(1, NbPrimitives);
  PrimLZ = dvector(1, NbPrimitives);
  Radius = dvector(1, NbPrimitives);  // can lead to a little memory misuse
  PrimOriginX = dvector(1, NbPrimitives);
  PrimOriginY = dvector(1, NbPrimitives);
  PrimOriginZ = dvector(1, NbPrimitives);
  PrimDC = (DirnCosn3D *)malloc((NbPrimitives + 1) * sizeof(DirnCosn3D));
  VolRef1 = ivector(1, NbPrimitives);
  VolRef2 = ivector(1, NbPrimitives);
  NbSurfSegX = ivector(1, NbPrimitives);
  NbSurfSegZ = ivector(1, NbPrimitives);
  NbWireSeg = ivector(1, NbPrimitives);  // little memory misuse
  InterfaceType = ivector(1, NbPrimitives);
  Epsilon1 = dvector(1, NbPrimitives);
  Epsilon2 = dvector(1, NbPrimitives);
  Lambda = dvector(1, NbPrimitives);
  ApplPot = dvector(1, NbPrimitives);
  ApplCh = dvector(1, NbPrimitives);
  PeriodicTypeX = ivector(1, NbPrimitives);
  PeriodicTypeY = ivector(1, NbPrimitives);
  PeriodicTypeZ = ivector(1, NbPrimitives);
  PeriodicInX = ivector(1, NbPrimitives);
  PeriodicInY = ivector(1, NbPrimitives);
  PeriodicInZ = ivector(1, NbPrimitives);
  XPeriod = dvector(1, NbPrimitives);
  YPeriod = dvector(1, NbPrimitives);
  ZPeriod = dvector(1, NbPrimitives);
  MirrorTypeX = ivector(1, NbPrimitives);
  MirrorTypeY = ivector(1, NbPrimitives);
  MirrorTypeZ = ivector(1, NbPrimitives);
  MirrorDistXFromOrigin = dvector(1, NbPrimitives);
  MirrorDistYFromOrigin = dvector(1, NbPrimitives);
  MirrorDistZFromOrigin = dvector(1, NbPrimitives);
  BndPlaneInXMin = ivector(1, NbPrimitives);
  BndPlaneInXMax = ivector(1, NbPrimitives);
  BndPlaneInYMin = ivector(1, NbPrimitives);
  BndPlaneInYMax = ivector(1, NbPrimitives);
  BndPlaneInZMin = ivector(1, NbPrimitives);
  BndPlaneInZMax = ivector(1, NbPrimitives);
  XBndPlaneInXMin = dvector(1, NbPrimitives);
  XBndPlaneInXMax = dvector(1, NbPrimitives);
  YBndPlaneInYMin = dvector(1, NbPrimitives);
  YBndPlaneInYMax = dvector(1, NbPrimitives);
  ZBndPlaneInZMin = dvector(1, NbPrimitives);
  ZBndPlaneInZMax = dvector(1, NbPrimitives);
  VBndPlaneInXMin = dvector(1, NbPrimitives);
  VBndPlaneInXMax = dvector(1, NbPrimitives);
  VBndPlaneInYMin = dvector(1, NbPrimitives);
  VBndPlaneInYMax = dvector(1, NbPrimitives);
  VBndPlaneInZMin = dvector(1, NbPrimitives);
  VBndPlaneInZMax = dvector(1, NbPrimitives);
  NbElmntsOnPrim = ivector(1, NbPrimitives);
  ElementBgn = ivector(1, NbPrimitives);
  ElementEnd = ivector(1, NbPrimitives);
  AvChDen = dvector(1, NbPrimitives);
  AvAsgndChDen = dvector(1, NbPrimitives);
 
  // Loop over the primitives - major loop
  int nvertex, volref1, volref2, volmax = 0;
#ifdef __cplusplus
  std::vector<double> xvert(MaxNbVertices, 0.);
  std::vector<double> yvert(MaxNbVertices, 0.);
  std::vector<double> zvert(MaxNbVertices, 0.);
#else
  double xvert[MaxNbVertices], yvert[MaxNbVertices], zvert[MaxNbVertices];
#endif
  double xnorm, ynorm, znorm;  // in case of wire , radius is read as xnorm
  for (int prim = 1; prim <= NbPrimitives; ++prim) {
#ifdef __cplusplus
    fstatus = neBEMGetPrimitive(prim, &nvertex, xvert.data(), yvert.data(),
                                zvert.data(), &xnorm, &ynorm, &znorm, &volref1,
                                &volref2);
#else
    fstatus = neBEMGetPrimitive(prim, &nvertex, xvert, yvert, zvert,  // arrays
                                &xnorm, &ynorm, &znorm, &volref1, &volref2);
#endif
    if (fstatus != 0) {
      neBEMMessage("neBEMReadGeometry - neBEMGetPrimitve");
      return -1;
    }
    if (volmax < volref1) {
      volmax = volref1;
    }  // maxm nb of volumes
    if (volmax < volref2) {
      volmax = volref2;
    }  // maxm nb of volumes
 
    if (nvertex > MaxNbVertices) {
      printf("Number of vertices for primitive %d exceeds %d!\n", prim,
             MaxNbVertices);
      printf("Returning to garfield ...\n");
      return (-1);
    }
 
    PrimType[prim] = nvertex;  // wire:2, triangle:3, rectangle:4
    NbVertices[prim] = nvertex;
    for (int vert = 0; vert < NbVertices[prim]; ++vert) {
      XVertex[prim][vert] = xvert[vert];
      YVertex[prim][vert] = yvert[vert];
      ZVertex[prim][vert] = zvert[vert];
    }
    if (PrimType[prim] == 2) {
      // wire
      XNorm[prim] = 0.0;  // modulus not 1 - an absurd trio!
      YNorm[prim] = 0.0;
      ZNorm[prim] = 0.0;
      Radius[prim] = xnorm;
    }
    if ((PrimType[prim] == 3) || (PrimType[prim] == 4)) {
      XNorm[prim] = xnorm;
      YNorm[prim] = ynorm;
      ZNorm[prim] = znorm;
      Radius[prim] = 0.0;  // absurd radius!
    }
    VolRef1[prim] = volref1;
    VolRef2[prim] = volref2;
 
    // feedback for user begins - suppress later
    if (OptPrintPrimaryDetails) {
      printf("neBEM:\tprimitive %d between volumes %d, %d has %d vertices\n",
             prim, volref1, volref2, nvertex);
      for (int ivertex = 0; ivertex < nvertex; ivertex++) {
        printf("\tnode %d (%g,%g,%g)\n", ivertex, xvert[ivertex],
               yvert[ivertex], zvert[ivertex]);
      }
      printf("\tnormal vector: (%g,%g,%g)\n", xnorm, ynorm, znorm);
    }
    // feedback for user ends - suppress later
 
    // Now look for the volume related information for this primitve
    // This is obtained from the specified volume references
    // volref1 refers to the volume itself
    // volref2 describes the volume in the direction of the +ve normal
    // Note that materials from 1 to 10 are conductors and
    //                                       from 11 to 20
    // are dielectrics
    int shape1, material1, boundarytype1;
    double epsilon1, potential1, charge1;
    if (volref1 == -1) {
      // Must be an error, since no device is made of vacuum
      neBEMMessage("neBEMReadGeometry - volref1 = -1!");
      return -1;
    } else {
      neBEMVolumeDescription(volref1, &shape1, &material1, &epsilon1,
                             &potential1, &charge1, &boundarytype1);
    }
    if (OptPrintVolumeDetails) {
      printf("\tvolref1: %d\n", volref1);
      printf("\t\tboundarytype1: %d, shape1: %d, material1: %d\n",
             boundarytype1, shape1, material1);
      printf("\t\tepsilon1: %lg, potential1: %lg, charge1: %lg\n", epsilon1,
             potential1, charge1);
    }
    // in the -ve normal direction - properties of the external volume
    int shape2, material2, boundarytype2;
    double epsilon2, potential2, charge2;
    if (volref2 == -1) {
      shape2 = 0;
      material2 = 11;
      epsilon2 = 1.0;
      potential2 = 0.0;
      charge2 = 0.0;
      boundarytype2 = 0;
    } else {
      neBEMVolumeDescription(volref2, &shape2, &material2, &epsilon2,
                             &potential2, &charge2, &boundarytype2);
    }
    if (OptPrintVolumeDetails) {
      printf("\tvolref2: %d\n", volref2);
      printf("\t\tboundarytype2: %d, shape2: %d, material2: %d\n",
             boundarytype2, shape2, material2);
      printf("\t\tepsilon2: %lg, potential2: %lg, charge2: %lg\n", epsilon2,
             potential2, charge2);
    }
 
    // Put default values to variables that depend on the interface type
    // Is there any risk involved in putting these defaults?
    // At present, they even seem necessary. For example, for floating
    // conductors or dielectric-dielectric interface, the formulation requires
    // that the RHS is zero (may be modified by the effect of known charges).
    Epsilon1[prim] = epsilon1;
    Epsilon2[prim] = epsilon2;  // 1: self, 2: external
    ApplPot[prim] = 0.0;
    Lambda[prim] = 0.0;
    ApplCh[prim] = 0.0;
 
    // BoundaryTypes:
    // --------------
    // Vacuum: 0
    // Conductor at specified potential: 1
    // Conductor with a specified charge: 2
    // Floating conductor (zero charge, perpendicular E): 3
    // Dielectric intrface (plastic-plastic) without "manual" charge: 4
    // Dielectric with surface charge (plastic-gas, typically): 5
    // Symmetry boundary, E parallel: 6 (may not be necessary)
    // Symmetry boundary, E perpendicular: 7 (may not be necessary)
 
    // InterfaceType:
    // --------------
    // To be skipped: 0
    // Conductor-dielectric: 1
    // Conductor with known charge: 2
    // Conductor at floating potential: 3
    // Dielectric-dielectric: 4
    // Dielectric with given surface charge: 5
    // Check dielectric-dielectric formulation in
    // (NumSolnOfBIEforMolES_JPBardhan.pdf):
    // Numerical solution of boundary-integral equations for molecular
    // electrostatics,
    // by Jaydeep P. Bardhan,
    // THE JOURNAL OF CHEMICAL PHYSICS 130, 094102 (2009)
 
    switch (boundarytype1) {  // the volume itself is volref1
      case 1:                 // conductor at specified potential
        if (boundarytype2 == 0 || boundarytype2 == 4) {
          // dielectric-conductor
          InterfaceType[prim] = 1;
          ApplPot[prim] = potential1;
        } else if (boundarytype2 == 1) {
          // conductor-conductor
          if (fabs(potential1 - potential2)  // same potential
              < 1e-6 * (1 + fabs(potential1) + fabs(potential2))) {
            printf("neBEMReadGeometry: identical potentials; skipped.\n");
            printf("Primitive skipped: #%d\n", prim);
            InterfaceType[prim] = 0;
          } else {
            // different potentials
            printf("neBEMReadGeometry: different potentials; rejected.\n");
            return -1;
          }
        } else {
          // conductor-unknown
          printf(
              "neBEMReadGeometry: conductor at given potential; rejected.\n");
          return -1;
        }
        break;
 
      case 2:  // conductor with a specified charge
        if ((boundarytype2 == 0) || (boundarytype2 == 4)) {
          // conductor-dielectric
          InterfaceType[prim] = 2;
          ApplCh[prim] = charge1;
        } else {
          printf("neBEMReadGeometry: charged conductor; rejected.\n");
          return -1;
        }
        break;
 
      case 3:  // floating conductor (zero charge, perpendicular E)
        if ((boundarytype2 == 0) || (boundarytype2 == 4)) {
          // conductor-dielectric
          InterfaceType[prim] = 3;
          if (!NbFloatingConductors)   // assuming only one floating conductor
            NbFloatingConductors = 1;  // in the system
        } else {
          printf("neBEMReadGeometry: floating conductor; rejected.\n");
          return -1;
        }
        break;
 
      case 4:  // dielectric interface (plastic-plastic) without "manual" charge
        if (boundarytype2 == 0) {
          // dielectric-vacuum
          // epsilon1 is self dielectric-constant
          // epsilon2 is towards positive normal
          InterfaceType[prim] = 4;
          Lambda[prim] = (epsilon1 - epsilon2) / (epsilon1 + epsilon2);
          // consistent with Bardhan's eqn 16 where (1 / (2*Lambda)) is used
        } else if (boundarytype2 == 1) {
          // dielectric-conductor
          InterfaceType[prim] = 1;  // conductor at known potential
          ApplPot[prim] = potential2;
        } else if (boundarytype2 == 2) {
          // dielectric-conductor
          InterfaceType[prim] = 2;  // conductor with known charge
          ApplCh[prim] = charge2;
        } else if (boundarytype2 == 3) {
          // dielectric-conductor
          InterfaceType[prim] = 3;  // conductor at floating potential
        } else if (boundarytype2 == 4) {
          // dielectric-dielectric
          if (fabs(epsilon1 - epsilon2) <
              1e-6 * (1 + fabs(epsilon1) + fabs(epsilon2))) {
            // identical dielectrica
            printf(
                "neBEMReadGeometry: between identical dielectrica; skipd.\n");
            printf("Primitive skipped: #%d\n", prim);
            InterfaceType[prim] = 0;
          } else {
            // distinctly different dielectrica
            // epsilon1 is self dielectric-constant
            // epsilon2 towards positive normal
            InterfaceType[prim] = 4;
            Lambda[prim] = (epsilon1 - epsilon2) / (epsilon1 + epsilon2);
            // consistent with Bardhan's paper (1 / Lambda)
          }
        } else if (boundarytype2 == 5) {
          // dielectric-dielectric with charge
          if (fabs(epsilon1 - epsilon2)  // identical dielectrica
              < 1e-6 * (1 + fabs(epsilon1) + fabs(epsilon2))) {
            printf(
                "neBEMReadGeometry: between identical dielectrica; skipped.\n");
            printf("Primitive skipped: #%d\n", prim);
            InterfaceType[prim] = 0;
          } else {
            // distinctly different dielectrica
            InterfaceType[prim] = 5;  // epsilon2 towards positive normal
            ApplCh[prim] = charge2;
            Lambda[prim] = (epsilon1 - epsilon2) / (epsilon1 + epsilon2);
          }
        }       // if-else if boundarytypes 0 and 4
        else {  // dielectric-unknown
          printf("neBEMReadGeometry: unknown dielectric; rejected.\n");
          return -1;
        }
        break;
 
      case 5:  // dielectric with surface charge (plastic-gas, typically)
        if (boundarytype2 == 0) {   // dielectric-vacuum
          InterfaceType[prim] = 5;  // epsilon2 is towards +ve normal
          ApplCh[prim] = charge1;
          Lambda[prim] = (epsilon1 - epsilon2) / (epsilon1 + epsilon2);
          // consistent with Bardhan's paper (1 / Lambda)
        } else if (boundarytype2 == 4) {  // dielectric-dielectric
          if (fabs(epsilon1 - epsilon2) <
              1e-6 * (1 + fabs(epsilon1) + fabs(epsilon2))) {
            // identical dielectrica
            printf(
                "neBEMReadGeometry: between identical dielectrica; skipd.\n");
            printf("Primitive skipped: #%d\n", prim);
            InterfaceType[prim] = 0;
          } else {
            // distinctly different dielectrica
            InterfaceType[prim] = 5;  // epsilon2 towards positive normal
            ApplCh[prim] = charge1;
            Lambda[prim] = (epsilon1 - epsilon2) / (epsilon1 + epsilon2);
            // consistent with Bardhan's paper (1 / Lambda)
          }
        }  // if-else if boundarytypes 0 and 4
        else {
          printf(
              "neBEMReadGeometry: charged dielectric adjacent to a conductor; "
              "rejected.\n");
          return -1;
        }
        break;
 
      case 6:  // symmetry boundary, E parallel
        if (boundarytype2 == 0) {
          InterfaceType[prim] = 6;
        } else {
          printf("neBEMReadGeometry: E-parallel symmetry; rejected.\n");
          return -1;
        }
        break;
 
      case 7:  // symmetry boundary, E perpendicular
        if (boundarytype2 == 0) {
          InterfaceType[prim] = 7;
        } else {
          printf("neBEMReadGeometry: E-perpendicular symmetry; rejected.\n");
          return -1;
        }
        break;
 
      default:
        printf("neBEMReadGeometry: Boundary type 1: %d\n", boundarytype1);
        printf("neBEMReadGeometry: Boundary type 2: %d\n", boundarytype2);
        printf("neBEMReadGeometry:          out of range ... exiting.\n");
        return -1;
    }  // switch boundarytype1 ends
 
    if (OptPrintVolumeDetails) {
      printf(
          "\tType: %d, ApplPot: %lg, Epsilon1: %lg, Epsilon2: %lg, Lambda: "
          "%lg, ApplCh: %lg\n",
          InterfaceType[prim], ApplPot[prim], Epsilon1[prim], Epsilon2[prim],
          Lambda[prim], ApplCh[prim]);
    }
 
    // Read the periodicities
    // Note that mirror has been taken care of separately below
    // ix: PeriodicTypeX (1 - simple, 2 - mirror, 3 - axial, 4 - rotation)
    // jx: PeriodicInX (Number of copies internal to neBEM)
    // sx: XPeriod
    // NOTE: A change in this part of the algorithm is likely to affect
    // src/Solve/neBEM.c (LHMatrix) and
    // src/Solve/ComputeProperties.c (PFAtPoint and WtFldPFAtPoint)
    {
      int ix, iy, iz;
      int jx, jy, jz;
      double sx, sy, sz;
      fstatus = neBEMGetPeriodicities(prim, &ix, &jx, &sx, &iy, &jy, &sy, &iz,
                                      &jz, &sz);
      if (fstatus != 0) {
        neBEMMessage("neBEMReadGeometry - neBEMGetPeriodicities");
        return -1;
      }
      if (jx < 0) jx = 0;
      if (jy < 0) jy = 0;
      if (jz < 0) jz = 0;
 
      PeriodicTypeX[prim] = ix;
      PeriodicTypeY[prim] = iy;
      PeriodicTypeZ[prim] = iz;
      if (0) {
        printf("For primitive: %d\n", prim);
        printf("\tPeriodicTypeX: %d, PeriodicTypeY: %d, PeriodicTypeZ: %d\n",
               ix, iy, iz);
        printf("\tPeriodicInX: %d, PeriodicInY: %d, PeriodicInZ: %d\n", jx, jy,
               jz);
        printf("\tXPeriod: %lg, YPeriod: %lg, ZPeriod: %lg\n", sx, sy, sz);
      }
      if (ix > 0) {
        // These checks need to be done separately. Otherwise, there is
        // a possibility that non-zero values of PeriodicIn* and *Period
        // are used throughout the code despite PeriodicType* is 0
        PeriodicInX[prim] = jx;
        XPeriod[prim] = sx;
      } else {
        PeriodicInX[prim] = 0;
        XPeriod[prim] = 0.0;
      }
      if (iy > 0) {
        PeriodicInY[prim] = jy;
        YPeriod[prim] = sy;
      } else {
        PeriodicInY[prim] = 0;
        YPeriod[prim] = 0.0;
      }
      if (iz > 0) {
        PeriodicInZ[prim] = jz;
        ZPeriod[prim] = sz;
      } else {
        PeriodicInZ[prim] = 0;
        ZPeriod[prim] = 0.0;
      }
    }  // read periodicity information
 
    // Read mirror information
    // Mirror can be in perpendicular to any of the three cartesian axes
    // There can be more than one mirror at the same time
    // MirrorType (1 - charge density -ve of original, equivalent to method of
    // images, 2 - charge density same as original)
    {
      int ix, iy, iz;
      int jx, jy, jz;  // not used at present
      double sx, sy, sz;
      fstatus =
          neBEMGetMirror(prim, &ix, &jx, &sx, &iy, &jy, &sy, &iz, &jz, &sz);
      if (fstatus != 0) {
        neBEMMessage("neBEMReadGeometry - neBEMGetMirror");
        return -1;
      }
      if (jx < 0) jx = 0;
      if (jy < 0) jy = 0;
      if (jz < 0) jz = 0;
 
      MirrorTypeX[prim] = ix;
      MirrorTypeY[prim] = iy;
      MirrorTypeZ[prim] = iz;
      if (0) {
        printf("For primitive: %d\n", prim);
        printf("\tMirrorTypeX: %d, MirrorTypeY: %d, MirrorTypeZ: %d\n", ix, iy,
               iz);
        printf("\tNOT USED ==> MirrorInX: %d, MirrorInY: %d, MirrorInZ: %d\n",
               jx, jy, jz);
        printf("\tMirrorDistX: %lg, MirrorDistY: %lg, MirrorDistZ: %lg\n", sx,
               sy, sz);
        getchar();
      }
      if (ix > 0) {
        // printf("neBEMReadGeometry: Mirror have been requested.\n");
        MirrorDistXFromOrigin[prim] = sx;  // assumed to pass through the origin
      } else {
        MirrorDistXFromOrigin[prim] = 0.0;  // pass through the origin
      }
      if (iy > 0) {
        // printf("neBEMReadGeometry: Mirror have been requested.\n");
        MirrorDistYFromOrigin[prim] = sy;
      } else {
        MirrorDistYFromOrigin[prim] = 0.0;
      }
      if (iz > 0) {
        // printf("neBEMReadGeometry: Mirror have been requested.\n");
        MirrorDistZFromOrigin[prim] = sz;
      } else {
        MirrorDistZFromOrigin[prim] = 0.0;
      }
    }  // read mirror information
 
    // Information on bounding planes
    // ixmin=0: lower x-plane does not exist
    // ixmin=1: lower x-plane does exist
    // cxmin: coordinate of lower x-plane
    // vxmin: potential of lower x-plane
    // Similar for ixmax, iymin, iymax, izmin, izmax
    int ixmin, ixmax, iymin, iymax, izmin, izmax;
    double cxmin, cxmax, cymin, cymax, czmin, czmax;
    double vxmin, vxmax, vymin, vymax, vzmin, vzmax;
    fstatus = neBEMGetBoundingPlanes(
        &ixmin, &cxmin, &vxmin, &ixmax, &cxmax, &vxmax, &iymin, &cymin, &vymin,
        &iymax, &cymax, &vymax, &izmin, &czmin, &vzmin, &izmax, &czmax, &vzmax);
    if (fstatus != 0) {
      neBEMMessage("neBEMReadGeometry - neBEMGetBoundingPlanes");
      return -1;
    }
    BndPlaneInXMin[prim] = ixmin;
    BndPlaneInXMax[prim] = ixmax;
    BndPlaneInYMin[prim] = iymin;
    BndPlaneInYMax[prim] = iymax;
    BndPlaneInZMin[prim] = izmin;
    BndPlaneInZMax[prim] = izmax;
    if (ixmin) {
      XBndPlaneInXMin[prim] = cxmin;
      VBndPlaneInXMin[prim] = vxmin;
    } else {
      XBndPlaneInXMin[prim] = 0.0;
      VBndPlaneInXMin[prim] = 0.0;
    }
    if (ixmax > 0) {
      XBndPlaneInXMax[prim] = cxmax;
      VBndPlaneInXMax[prim] = vxmax;
    } else {
      XBndPlaneInXMax[prim] = 0.0;
      VBndPlaneInXMax[prim] = 0.0;
    }
    if (iymin > 0) {
      YBndPlaneInYMin[prim] = cymin;
      VBndPlaneInYMin[prim] = vymin;
    } else {
      YBndPlaneInYMin[prim] = 0.0;
      VBndPlaneInYMin[prim] = 0.0;
    }
    if (iymax > 0) {
      YBndPlaneInYMax[prim] = cymax;
      VBndPlaneInYMax[prim] = vymax;
    } else {
      YBndPlaneInYMax[prim] = 0.0;
      VBndPlaneInYMax[prim] = 0.0;
    }
    if (izmin > 0) {
      ZBndPlaneInZMin[prim] = czmin;
      VBndPlaneInZMin[prim] = vzmin;
    } else {
      ZBndPlaneInZMin[prim] = 0.0;
      VBndPlaneInZMin[prim] = 0.0;
    }
    if (izmax > 0) {
      ZBndPlaneInZMax[prim] = czmax;
      VBndPlaneInZMax[prim] = vzmax;
    } else {
      ZBndPlaneInZMax[prim] = 0.0;
      VBndPlaneInZMax[prim] = 0.0;
    }
  }  // loop over the primitives - major loop
 
  VolMax = volmax;              // Maximum nb of volumes in the problem
  volRef = ivector(0, VolMax);  // variables to store volume related infomration
  volShape = ivector(0, VolMax);
  volMaterial = ivector(0, VolMax);
  volEpsilon = dvector(0, VolMax);
  volPotential = dvector(0, VolMax);
  volCharge = dvector(0, VolMax);
  volBoundaryType = ivector(0, VolMax);
  for (int volref = 0; volref <= VolMax; ++volref) {
    neBEMVolumeDescription(volref, &volShape[volref], &volMaterial[volref],
                           &volEpsilon[volref], &volPotential[volref],
                           &volCharge[volref], &volBoundaryType[volref]);
    if (dbgFn) {
      printf("volref: %d\n", volref);
      printf("shape: %d,  material: %d\n", volShape[volref],
             volMaterial[volref]);
      printf("eps: %lg,  pot: %lg\n", volEpsilon[volref], volPotential[volref]);
      printf("q: %lg,  type: %d\n", volCharge[volref], volBoundaryType[volref]);
    }
  }
 
  // Ignore unnecessary primitives from the final count
  // Ideally, all the removal conditions for a primitive should be checked in
  // one loop and the list should be updated in one single go.
  {
    for (int prim = 1; prim <= NbPrimitives; ++prim) {
      OrgnlToEffPrim[prim][0] = prim;
      OrgnlToEffPrim[prim][1] = prim;
      OrgnlToEffPrim[prim][2] = prim;
    }
 
    {  // Skip primitive
      // Remove skipped primitives having InterfaceType == 0.
      // Also remove primitives having too small dimensions.
      int NbSkipped = 0, effprim;
      double DVertex[4], minDVertex = 0.0;  // maximum number of vertices is 4
      for (int prim = 1; prim <= NbPrimitives; ++prim) {
        effprim = prim - NbSkipped;
 
        // Check dimensions of the primitive
        for (int vert = 0; vert < NbVertices[prim] - 1; ++vert) {
          DVertex[vert] =
              sqrt(((XVertex[prim][vert + 1] - XVertex[prim][vert]) *
                    (XVertex[prim][vert + 1] - XVertex[prim][vert])) +
                   ((YVertex[prim][vert + 1] - YVertex[prim][vert]) *
                    (YVertex[prim][vert + 1] - YVertex[prim][vert])) +
                   ((ZVertex[prim][vert + 1] - ZVertex[prim][vert]) *
                    (ZVertex[prim][vert + 1] - ZVertex[prim][vert])));
          if (vert == 0)
            minDVertex = DVertex[vert];
          else {
            if (DVertex[vert] < minDVertex) minDVertex = DVertex[vert];
          }
        }
 
        if ((InterfaceType[prim]) && (minDVertex > MINDIST)) {
          OrgnlToEffPrim[prim][1] = effprim;
          OrgnlToEffPrim[prim][2] = effprim;
          PrimType[effprim] = PrimType[prim];
          NbVertices[effprim] = NbVertices[prim];
          for (int vert = 0; vert < NbVertices[effprim]; ++vert) {
            XVertex[effprim][vert] = XVertex[prim][vert];
            YVertex[effprim][vert] = YVertex[prim][vert];
            ZVertex[effprim][vert] = ZVertex[prim][vert];
          }
          if (PrimType[effprim] == 2)  // wire
          {
            XNorm[effprim] = 0.0;  // modulus not 1 - an absurd trio!
            YNorm[effprim] = 0.0;
            ZNorm[effprim] = 0.0;
            Radius[effprim] = Radius[prim];
          }
          if ((PrimType[effprim] == 3) || (PrimType[effprim] == 4)) {
            XNorm[effprim] = XNorm[prim];
            YNorm[effprim] = YNorm[prim];
            ZNorm[effprim] = ZNorm[prim];
            Radius[effprim] = 0.0;  // absurd radius!
          }
          VolRef1[effprim] = VolRef1[prim];
          VolRef2[effprim] = VolRef2[prim];
 
          InterfaceType[effprim] = InterfaceType[prim];
          Epsilon1[effprim] = Epsilon1[prim];
          Epsilon2[effprim] = Epsilon2[prim];
          Lambda[effprim] = Lambda[prim];
          ApplPot[effprim] = ApplPot[prim];
          ApplCh[effprim] = ApplCh[prim];
          PeriodicTypeX[effprim] = PeriodicTypeX[prim];
          PeriodicTypeY[effprim] = PeriodicTypeY[prim];
          PeriodicTypeZ[effprim] = PeriodicTypeZ[prim];
          PeriodicInX[effprim] = PeriodicInX[prim];
          PeriodicInY[effprim] = PeriodicInY[prim];
          PeriodicInZ[effprim] = PeriodicInZ[prim];
          XPeriod[effprim] = XPeriod[prim];
          YPeriod[effprim] = YPeriod[prim];
          ZPeriod[effprim] = ZPeriod[prim];
          MirrorTypeX[effprim] = MirrorTypeX[prim];
          MirrorTypeY[effprim] = MirrorTypeY[prim];
          MirrorTypeZ[effprim] = MirrorTypeZ[prim];
          MirrorDistXFromOrigin[effprim] = MirrorDistXFromOrigin[prim];
          MirrorDistYFromOrigin[effprim] = MirrorDistYFromOrigin[prim];
          MirrorDistZFromOrigin[effprim] = MirrorDistZFromOrigin[prim];
          BndPlaneInXMin[effprim] = BndPlaneInXMin[prim];
          BndPlaneInXMax[effprim] = BndPlaneInXMax[prim];
          BndPlaneInYMin[effprim] = BndPlaneInYMin[prim];
          BndPlaneInYMax[effprim] = BndPlaneInYMax[prim];
          BndPlaneInZMin[effprim] = BndPlaneInZMin[prim];
          BndPlaneInZMax[effprim] = BndPlaneInZMax[prim];
          XBndPlaneInXMin[effprim] = XBndPlaneInXMin[prim];
          XBndPlaneInXMax[effprim] = XBndPlaneInXMax[prim];
          YBndPlaneInYMin[effprim] = YBndPlaneInYMin[prim];
          YBndPlaneInYMax[effprim] = YBndPlaneInYMax[prim];
          ZBndPlaneInZMin[effprim] = ZBndPlaneInZMin[prim];
          ZBndPlaneInZMax[effprim] = ZBndPlaneInZMax[prim];
          VBndPlaneInXMin[effprim] = VBndPlaneInXMin[prim];
          VBndPlaneInXMax[effprim] = VBndPlaneInXMax[prim];
          VBndPlaneInYMin[effprim] = VBndPlaneInYMin[prim];
          VBndPlaneInYMax[effprim] = VBndPlaneInYMax[prim];
          VBndPlaneInZMin[effprim] = VBndPlaneInZMin[prim];
          VBndPlaneInZMax[effprim] = VBndPlaneInZMax[prim];
        }  // InterfaceType
        else {
          OrgnlToEffPrim[prim][1] = 0;  // removed from the list
          OrgnlToEffPrim[prim][2] = 0;
          ++NbSkipped;
          if (DebugLevel == 101) {
            printf("Skipped primitive %d, InterfaceType: %d, minDVertex: %lg\n",
                   prim, InterfaceType[prim], minDVertex);
          }
        }
      }  // loop over primitives to remove the skipped primitives
      NbPrimitives -= NbSkipped;
      printf("Number of primitives skipped: %d, Effective NbPrimitives: %d\n",
             NbSkipped, NbPrimitives);
    }  // Skip primitives
 
    if (OptRmPrim) {
      int NbRmPrims;
      FILE *rmprimFile = fopen("neBEMInp/neBEMRmPrim.inp", "r");
      if (rmprimFile == NULL) {
        printf("neBEMRmPrim.inp absent ... assuming defaults ...\n");
        NbRmPrims = 0;
      } else {
        fscanf(rmprimFile, "NbRmPrims: %d\n", &NbRmPrims);
        if (NbRmPrims) {
          int tint;
#ifdef __cplusplus
          std::vector<double> rmXNorm(NbRmPrims + 1, 0.);
          std::vector<double> rmYNorm(NbRmPrims + 1, 0.);
          std::vector<double> rmZNorm(NbRmPrims + 1, 0.);
          std::vector<double> rmXVert(NbRmPrims + 1, 0.);
          std::vector<double> rmYVert(NbRmPrims + 1, 0.);
          std::vector<double> rmZVert(NbRmPrims + 1, 0.);
#else
          double rmXNorm[NbRmPrims + 1], rmYNorm[NbRmPrims + 1];
          double rmZNorm[NbRmPrims + 1];
          double rmXVert[NbRmPrims + 1], rmYVert[NbRmPrims + 1];
          double rmZVert[NbRmPrims + 1];
#endif
          for (int rmprim = 1; rmprim <= NbRmPrims; ++rmprim) {
            fscanf(rmprimFile, "Prim: %d\n", &tint);
            fscanf(rmprimFile, "rmXNorm: %le\n", &rmXNorm[rmprim]);
            fscanf(rmprimFile, "rmYNorm: %le\n", &rmYNorm[rmprim]);
            fscanf(rmprimFile, "rmZNorm: %le\n", &rmZNorm[rmprim]);
            fscanf(rmprimFile, "rmXVert: %le\n", &rmXVert[rmprim]);
            fscanf(rmprimFile, "rmYVert: %le\n", &rmYVert[rmprim]);
            fscanf(rmprimFile, "rmZVert: %le\n", &rmZVert[rmprim]);
            printf(
                "rmprim: %d, rmXNorm: %lg, rmYNorm: %lg, rmZNorm: %lg, "
                "rmXVert: %lg, rmYVert: %lg, rmZVert: %lg\n",
                rmprim, rmXNorm[rmprim], rmYNorm[rmprim], rmZNorm[rmprim],
                rmXVert[rmprim], rmYVert[rmprim], rmZVert[rmprim]);
          }
#ifdef __cplusplus
          std::vector<int> remove(NbPrimitives + 1, 0);
#else
          int remove[NbPrimitives + 1];
#endif
          // Check updated prim list
          for (int prim = 1; prim <= NbPrimitives; ++prim) {
            remove[prim] = 0;
            if (dbgFn) {
              printf("\n\nprim: %d, XVertex: %lg, YVertex: %lg, ZVertex: %lg\n",
                     prim, XVertex[prim][0], YVertex[prim][0],
                     ZVertex[prim][0]);
              printf("XNorm: %lg, YNorm: %lg, ZNorm: %lg\n", XNorm[prim],
                     YNorm[prim], ZNorm[prim]);
            }
 
            for (int rmprim = 1; rmprim <= NbRmPrims; ++rmprim) {
              if (dbgFn) {
                printf(
                    "rmprim: %d, rmXVertex: %lg, rmYVertex: %lg, rmZVertex: "
                    "%lg\n",
                    rmprim, rmXVert[rmprim], rmYVert[rmprim], rmZVert[rmprim]);
                printf("rmXNorm: %lg, rmYNorm: %lg, rmZNorm: %lg\n",
                       rmXNorm[rmprim], rmYNorm[rmprim], rmZNorm[rmprim]);
              }
 
              // check the normal
              if ((fabs(fabs(XNorm[prim]) - fabs(rmXNorm[rmprim])) <=
                   MINDIST) &&
                  (fabs(fabs(YNorm[prim]) - fabs(rmYNorm[rmprim])) <=
                   MINDIST) &&
                  (fabs(fabs(ZNorm[prim]) - fabs(rmZNorm[rmprim])) <=
                   MINDIST)) {  // prim and rmprim are parallel
                // coplanarity check to be implemented later.
                // For the time-being, we will assume that the planes to be
                // removed have their normals parallel to a given axis. So, we
                // only check that and remove the primitive if the distace along
                // that axis match. Possible pitfall => the primitives may be
                // coplanar but non-overlapping!
                if (fabs(fabs(XNorm[prim]) - 1.0) <= 1.0e-12) {
                  // primitive || to YZ
                  if (fabs(XVertex[prim][0] - rmXVert[rmprim]) <= MINDIST) {
                    remove[prim] = 1;
                  }
                }
                if (fabs(fabs(YNorm[prim]) - 1.0) <= 1.0e-12) {
                  // primitive || to XZ
                  if (fabs(YVertex[prim][0] - rmYVert[rmprim]) <= MINDIST) {
                    remove[prim] = 1;
                  }
                }
                if (fabs(fabs(ZNorm[prim]) - 1.0) <= 1.0e-12) {
                  // primitive || to XY
                  if (fabs(ZVertex[prim][0] - rmZVert[rmprim]) <= MINDIST) {
                    remove[prim] = 1;
                  }
                }
              }  // case where prim and rmprim are parallel
              if (dbgFn) {
                printf("prim: %d, rmprim: %d, remove: %d\n", prim, rmprim,
                       remove[prim]);
              }
              if (remove[prim] == 1)
                break;  // once removed, no point checking others
            }           // for rmprim - loop over all removal specification
 
          }  // for prim loop over all primitives
 
          int NbRemoved = 0;
          char RmPrimFile[256];
          strcpy(RmPrimFile, NativePrimDir);
          strcat(RmPrimFile, "/RmPrims.info");
          FILE *fprrm = fopen(RmPrimFile, "w");
          if (fprrm == NULL) {
            printf(
                "error opening RmPrims.info file in write mode ... "
                "returning\n");
 
            fclose(rmprimFile);
            return (-1);
          }
          // Note that some of the original primitives have already been removed
          // based on interface and dimension considerations
          int orgnlNb = 0;
          for (int prim = 1; prim <= NbPrimitives; ++prim) {
            // identify primitive number in the original list
            for (int orgnlprim = 1; orgnlprim <= OrgnlNbPrimitives;
                 ++orgnlprim) {
              if (OrgnlToEffPrim[orgnlprim][1] ==
                  prim)  // number updated for intrfc
              {
                orgnlNb = orgnlprim;
                break;
              }
            }  // loop for finding out its position in the original list
 
            if (remove[prim] == 1) {
              ++NbRemoved;
              OrgnlToEffPrim[orgnlNb][2] = 0;
              fprintf(fprrm, "NbRemoved: %d, Removed primitive: %d\n",
                      NbRemoved, prim);
              fprintf(fprrm, "PrimType: %d\n", PrimType[prim]);
              fprintf(fprrm, "NbVertices: %d\n", NbVertices[prim]);
              for (int vert = 0; vert < NbVertices[prim]; ++vert) {
                fprintf(fprrm, "Vertx %d: %lg, %lg, %lg\n", vert,
                        XVertex[prim][vert], YVertex[prim][vert],
                        ZVertex[prim][vert]);
              }
              fprintf(fprrm, "Normals: %lg, %lg, %lg\n", XNorm[prim],
                      YNorm[prim], ZNorm[prim]);
              continue;
            }       // if remove
            else {  // keep this one in the updated list of primitives
              int effprim = prim - NbRemoved;
 
              OrgnlToEffPrim[orgnlNb][2] = effprim;
              PrimType[effprim] = PrimType[prim];
              NbVertices[effprim] = NbVertices[prim];
              for (int vert = 0; vert < NbVertices[effprim]; ++vert) {
                XVertex[effprim][vert] = XVertex[prim][vert];
                YVertex[effprim][vert] = YVertex[prim][vert];
                ZVertex[effprim][vert] = ZVertex[prim][vert];
              }
              if (PrimType[effprim] == 2) {
                // wire
                XNorm[effprim] = 0.0;  // modulus not 1 - an absurd trio!
                YNorm[effprim] = 0.0;
                ZNorm[effprim] = 0.0;
                Radius[effprim] = Radius[prim];
              }
              if ((PrimType[effprim] == 3) || (PrimType[effprim] == 4)) {
                XNorm[effprim] = XNorm[prim];
                YNorm[effprim] = YNorm[prim];
                ZNorm[effprim] = ZNorm[prim];
                Radius[effprim] = 0.0;  // absurd radius!
              }
              VolRef1[effprim] = VolRef1[prim];
              VolRef2[effprim] = VolRef2[prim];
 
              InterfaceType[effprim] = InterfaceType[prim];
              Epsilon1[effprim] = Epsilon1[prim];
              Epsilon2[effprim] = Epsilon2[prim];
              Lambda[effprim] = Lambda[prim];
              ApplPot[effprim] = ApplPot[prim];
              ApplCh[effprim] = ApplCh[prim];
              PeriodicTypeX[effprim] = PeriodicTypeX[prim];
              PeriodicTypeY[effprim] = PeriodicTypeY[prim];
              PeriodicTypeZ[effprim] = PeriodicTypeZ[prim];
              PeriodicInX[effprim] = PeriodicInX[prim];
              PeriodicInY[effprim] = PeriodicInY[prim];
              PeriodicInZ[effprim] = PeriodicInZ[prim];
              XPeriod[effprim] = XPeriod[prim];
              YPeriod[effprim] = YPeriod[prim];
              ZPeriod[effprim] = ZPeriod[prim];
              MirrorTypeX[effprim] = MirrorTypeX[prim];
              MirrorTypeY[effprim] = MirrorTypeY[prim];
              MirrorTypeZ[effprim] = MirrorTypeZ[prim];
              MirrorDistXFromOrigin[effprim] = MirrorDistXFromOrigin[prim];
              MirrorDistYFromOrigin[effprim] = MirrorDistYFromOrigin[prim];
              MirrorDistZFromOrigin[effprim] = MirrorDistZFromOrigin[prim];
              BndPlaneInXMin[effprim] = BndPlaneInXMin[prim];
              BndPlaneInXMax[effprim] = BndPlaneInXMax[prim];
              BndPlaneInYMin[effprim] = BndPlaneInYMin[prim];
              BndPlaneInYMax[effprim] = BndPlaneInYMax[prim];
              BndPlaneInZMin[effprim] = BndPlaneInZMin[prim];
              BndPlaneInZMax[effprim] = BndPlaneInZMax[prim];
              XBndPlaneInXMin[effprim] = XBndPlaneInXMin[prim];
              XBndPlaneInXMax[effprim] = XBndPlaneInXMax[prim];
              YBndPlaneInYMin[effprim] = YBndPlaneInYMin[prim];
              YBndPlaneInYMax[effprim] = YBndPlaneInYMax[prim];
              ZBndPlaneInZMin[effprim] = ZBndPlaneInZMin[prim];
              ZBndPlaneInZMax[effprim] = ZBndPlaneInZMax[prim];
              VBndPlaneInXMin[effprim] = VBndPlaneInXMin[prim];
              VBndPlaneInXMax[effprim] = VBndPlaneInXMax[prim];
              VBndPlaneInYMin[effprim] = VBndPlaneInYMin[prim];
              VBndPlaneInYMax[effprim] = VBndPlaneInYMax[prim];
              VBndPlaneInZMin[effprim] = VBndPlaneInZMin[prim];
              VBndPlaneInZMax[effprim] = VBndPlaneInZMax[prim];
            }  // else remove == 0
          }    // loop over primitives to remove the primitives tagged to be
               // removed
          fclose(fprrm);
 
          NbPrimitives -= NbRemoved;
          printf(
              "Number of primitives removed: %d, Effective NbPrimitives: %d\n",
              NbRemoved, NbPrimitives);
          fflush(stdout);
        }  // if NbRmPrims true, implying primitives need to be removed
        fclose(rmprimFile);
      }  // if the rmprimFile is not NULL, prepare to remove primitives
    }    // if OptRmPrim: remove primitives as desired by the user
 
    // Information about primitives which are being ignored
    char IgnorePrimFile[256];
    strcpy(IgnorePrimFile, NativePrimDir);
    strcat(IgnorePrimFile, "/IgnorePrims.info");
    FILE *fignore = fopen(IgnorePrimFile, "w");
    if (fignore == NULL) {
      printf(
          "error opening IgnorePrims.info file in write mode ... returning\n");
      return (-1);
    }
 
    for (int prim = 1; prim <= OrgnlNbPrimitives; ++prim) {
      fprintf(fignore, "%d %d %d\n", OrgnlToEffPrim[prim][0],
              OrgnlToEffPrim[prim][1], OrgnlToEffPrim[prim][2]);
    }
 
    fclose(fignore);
  }  // Ignore unnecessary primitives from the final count
 
  // Reduced-Order Modelling information
  printf("ROM: switch to primitive representation after %d repetitions.\n",
         PrimAfter);
 
  // Store model data in native neBEM format
  char NativeInFile[256];
 
  strcpy(NativeInFile, NativeOutDir);
  strcat(NativeInFile, "/neBEMNative.inp");
  FILE *fNativeInFile = fopen(NativeInFile, "w");
  fprintf(fNativeInFile, "#====>Input directory\n");
  fprintf(fNativeInFile, "%s\n", NativePrimDir);
  fprintf(fNativeInFile, "#====>No. of primitives:\n");
  fprintf(fNativeInFile, "%d\n", NbPrimitives);
  fprintf(fNativeInFile, "#====>No. of volumes:\n");
  fprintf(fNativeInFile, "%d\n", VolMax);
  for (int prim = 1; prim <= NbPrimitives; ++prim) {
    char strPrimNb[11];
    snprintf(strPrimNb, 11, "%d", prim);
    char NativePrimFile[256];
    strcpy(NativePrimFile, "Primitive");
    strcat(NativePrimFile, strPrimNb);
    strcat(NativePrimFile, ".inp");
 
    fprintf(fNativeInFile, "#Input filename for primitive:\n");
    fprintf(fNativeInFile, "%s\n", NativePrimFile);
 
    strcpy(NativePrimFile, NativePrimDir);
    strcat(NativePrimFile, "/Primitive");
    strcat(NativePrimFile, strPrimNb);
    strcat(NativePrimFile, ".inp");
 
    FILE *fNativePrim = fopen(NativePrimFile, "w");
 
    fprintf(fNativePrim, "#Number of vertices:\n");
    fprintf(fNativePrim, "%d\n", NbVertices[prim]);
    for (int vertex = 0; vertex < NbVertices[prim]; ++vertex) {
      fprintf(fNativePrim, "#Vertex %d:\n", vertex);
      fprintf(fNativePrim, "%lg %lg %lg\n", XVertex[prim][vertex],
              YVertex[prim][vertex], ZVertex[prim][vertex]);
    }  // for vertex
    fprintf(fNativePrim, "#Outward normal:\n");
    fprintf(fNativePrim, "%lg %lg %lg\n", XNorm[prim], YNorm[prim],
            ZNorm[prim]);
    fprintf(fNativePrim, "#Volume references:\n");
    fprintf(fNativePrim, "%d %d\n", VolRef1[prim], VolRef2[prim]);
    fprintf(fNativePrim, "#Maximum number of segments:\n");
    fprintf(fNativePrim, "%d %d\n", 0, 0);  // use the trio target, min, max
    fprintf(fNativePrim, "#Information on X periodicity:\n");
    fprintf(fNativePrim, "%d %d %lg\n", PeriodicTypeX[prim], PeriodicInX[prim],
            XPeriod[prim]);
    fprintf(fNativePrim, "#Information on Y periodicity:\n");
    fprintf(fNativePrim, "%d %d %lg\n", PeriodicTypeY[prim], PeriodicInY[prim],
            YPeriod[prim]);
    fprintf(fNativePrim, "#Information on Z periodicity:\n");
    fprintf(fNativePrim, "%d %d %lg\n", PeriodicTypeZ[prim], PeriodicInZ[prim],
            ZPeriod[prim]);
 
    fclose(fNativePrim);
  }  // for prim
 
  fprintf(fNativeInFile, "#====>No. of boundary conditions per element:\n");
  fprintf(fNativeInFile, "%d\n", 1);  // CHECK!!!
  fprintf(fNativeInFile, "#====>Device output directory name:\n");
  fprintf(fNativeInFile, "NativeResults\n");
  fprintf(fNativeInFile, "#====>Map input file:\n");
  fprintf(fNativeInFile, "MapFile.inp\n");
  fclose(fNativeInFile);
 
  for (int volume = 0; volume <= VolMax; ++volume) {
    // Note that materials from 1 to 10 are conductors and
    //                     from 11 to 20 are dielectrics
    int shape, material, boundarytype;
    double epsilon, potential, charge;
    neBEMVolumeDescription(volume, &shape, &material, &epsilon, &potential,
                           &charge, &boundarytype);
 
    char strVolNb[11];
    snprintf(strVolNb, 11, "%d", volume);
    char NativeVolFile[256];
    strcpy(NativeVolFile, NativePrimDir);
    strcat(NativeVolFile, "/Volume");
    strcat(NativeVolFile, strVolNb);
    strcat(NativeVolFile, ".inp");
 
    FILE *fNativeVol = fopen(NativeVolFile, "w");
 
    fprintf(fNativeVol, "#Shape of the volume:\n");
    fprintf(fNativeVol, "%d\n", shape);
    fprintf(fNativeVol, "#Material of the volume:\n");
    fprintf(fNativeVol, "%d\n", material);
    fprintf(fNativeVol, "#Relative permittivity:\n");
    fprintf(fNativeVol, "%lg\n", epsilon);
    fprintf(fNativeVol, "#Potential:\n");
    fprintf(fNativeVol, "%lg\n", potential);
    fprintf(fNativeVol, "#Charge:\n");
    fprintf(fNativeVol, "%lg\n", charge);
    fprintf(fNativeVol, "#Boundary type:\n");
    fprintf(fNativeVol, "%d\n", boundarytype);
 
    fclose(fNativeVol);
  }
 
  // create a dummy map file
  {
    char NativeMapFile[256];
 
    strcpy(NativeMapFile, NativePrimDir);
    strcat(NativeMapFile, "/MapFile.inp");
 
    FILE *fNativeMap = fopen(NativeMapFile, "w");
 
    fprintf(fNativeMap, "#Number of maps\n");
    fprintf(fNativeMap, "%d\n", 0);
    fprintf(fNativeMap, "#Number of lines\n");
    fprintf(fNativeMap, "%d\n", 0);
    fprintf(fNativeMap, "#Number of points\n");
    fprintf(fNativeMap, "%d\n", 0);
 
    fclose(fNativeMap);
  }
 
  // Store primitive related data in a file for a new model, if opted for
  if (NewModel && OptStorePrimitives) {
    if (OptFormattedFile) {
      fstatus = WritePrimitives();
      if (fstatus) {
        neBEMMessage("neBEMReadGeometry - problem writing Primtives.\n");
        return -1;
      }
    }  // formatted file
 
    if (OptUnformattedFile) {
      neBEMMessage(
          "neBEMReadGeometry - unformatted write not inplemented yet.\n");
      return -1;
    }  // unformatted file
  }    // store primitives
 
  printf("neBEMReadGeometry: Geometry read!\n");
 
  neBEMState = 3;  // info about primitives read in after initialization and Nbs
 
  stopClock = clock();
  neBEMTimeElapsed(startClock, stopClock);
  printf("to read geometry\n");
 
  return (0);  // Success!
}  // neBEMReadGeometry ends

neBEMSetDefaults()

INTFACEGLOBAL int neBEMSetDefaults

(

void

)

neBEMSolve()

INTFACEGLOBAL int neBEMSolve

(

void

)

Definition at line 1811 of file neBEMInterface.c.

                     {
  int dbgFn = 0;
 
  clock_t startSolveClock = clock();
 
  if (TimeStep < 1) {
    neBEMMessage("neBEMSolve - TimeStep cannot be less than one!;\n");
    neBEMMessage("             Please put TimeStep = 1 for static problems.\n");
  }
 
  if ((neBEMState == 5) || (neBEMState == 8)) {
    if (neBEMState == 8)  // neBEMState 8 must have inverted flag on
    {  // so it must be a case looking for solution with a NewBC
      if (NewBC == 0) {
        neBEMMessage("neBEMSolve - NewBC zero for neBEMState = 8!");
        neBEMMessage("           - Nothing to be done ... returning.");
        return -1;
      }
    }
 
    if (NewModel) {  // effectively, NewMesh = NewBC = NewPP = 1;
      int fstatus = ComputeSolution();
      if (fstatus != 0) {
        neBEMMessage("neBEMSolve - NewModel");
        return -1;
      }
    } else {  // NewModel == 0
      if (NewMesh) {
        // effectively, NewBC = NewPP = 1;
        int fstatus = ComputeSolution();
        if (fstatus != 0) {
          neBEMMessage("neBEMSolve - NewMesh");
          return -1;
        }
      } else {        // NewModel == NewMesh == 0
        if (NewBC) {  // effectively, NewPP = 1;
          int fstatus = ComputeSolution();
          if (fstatus != 0) {
            neBEMMessage("neBEMSolve - Failure computing new solution");
            return -1;
          }
        } else {  // NewBC == 0
          if (NewPP) {
            int fstatus = ReadSolution();
            if (fstatus != 0) {
              neBEMMessage("neBEMSolve - Failure reading solution");
              return (-1);
            }
          } else {  // NewPP == 0
            printf("neBEMSolve: Nothing to do ... returning ...\n");
            return (-1);
          }  // NewPP == 0
        }    // NewBC == 0
      }      // NewModel == NewDiscretization == 0
    }        // NewModel == 0
 
    neBEMState = 9;
  } else {
    printf("neBEMSolve: neBEMSolve can be called only in state 5 / 8 ...\n");
    printf("returning ...\n");
    return (-1);
  }
 
  if (FailureCntr) {
    printf(
        "neBEMSolve: Approximations were made while computing the influence "
        "coefficients.\n");
    printf("            Please check the \"%s/Isles.log\" file.\n", PPOutDir);
  }
 
  clock_t stopSolveClock = clock();
  neBEMTimeElapsed(startSolveClock, stopSolveClock);
  printf("to complete solve\n");
 
  // Prepare voxelized data that will be exported to Garfield++
  if (OptVoxel) {
    clock_t startVoxelClock = clock();
 
    int fstatus = VoxelFPR();
    if (fstatus != 0) {
      neBEMMessage("neBEMSolve - Failure computing VoxelFPR");
      return -1;
    }
 
    clock_t stopVoxelClock = clock();
    neBEMTimeElapsed(startVoxelClock, stopVoxelClock);
    printf("to compute VoxelFPR\n");
  }
 
  // Prepare 3dMap data that will be exported to Garfield++
  if (OptMap) {
    clock_t startMapClock = clock();
 
    int fstatus = MapFPR();
    if (fstatus != 0) {
      neBEMMessage("neBEMSolve - Failure computing MapFPR");
      return -1;
    }
 
    clock_t stopMapClock = clock();
    neBEMTimeElapsed(startMapClock, stopMapClock);
    printf("to compute MapFPR\n");
  }
 
  // allocate memory for potential and field components within FAST volume mesh
  // and compute / read FastVol data
  // Similar allocation, computation and reading may be necessary for the KnCh
  // effects.
  // The other approach could be to create fast volumes that always have both
  // the influences (elements + known charges) added together. This approach
  // seems more managable now and is being followed.
  // Thus, if we want to inspect the effects of elements and known charges
  // separately, we will have to generate one fast volume with OptKnCh = 0,
  // and another with OptKnCh = 1. Subtraction of these two fast volumes will
  // provide us with the effect of KnCh.
  if (OptFastVol) {
    int MaxXCells = BlkNbXCells[1];
    int MaxYCells = BlkNbYCells[1];
    int MaxZCells = BlkNbZCells[1];
    clock_t startFastClock = clock();
 
    // find maximum number of Xcells etc in all the blocks
    // simplifies memory allocation using nrutils but hogs memory!
    for (int block = 1; block <= FastVol.NbBlocks; ++block) {
      if (block == 1) {
        MaxXCells = BlkNbXCells[1];
        MaxYCells = BlkNbYCells[1];
        MaxZCells = BlkNbZCells[1];
      } else {
        if (MaxXCells < BlkNbXCells[block]) MaxXCells = BlkNbXCells[block];
        if (MaxYCells < BlkNbYCells[block]) MaxYCells = BlkNbYCells[block];
        if (MaxZCells < BlkNbZCells[block]) MaxZCells = BlkNbZCells[block];
      }
    }  // loop block for finding maxm cells among all the blocks
 
    if (dbgFn) {
      printf("OptFastVol: %d\n", OptFastVol);
      printf("NbPtSkip: %d\n", NbPtSkip);
      printf("OptStaggerFastVol: %d\n", OptStaggerFastVol);
      printf("NbStgPtSkip: %d\n", NbStgPtSkip);
      printf("OptReadFastPF: %d\n", OptReadFastPF);
      printf("LX: %le\n", FastVol.LX);
      printf("LY: %le\n", FastVol.LY);
      printf("LZ: %le\n", FastVol.LZ);
      printf("CornerX: %le\n", FastVol.CrnrX);
      printf("CornerY: %le\n", FastVol.CrnrY);
      printf("CornerZ: %le\n", FastVol.CrnrZ);
      printf("YStagger: %le\n", FastVol.YStagger);
      printf("NbOfBlocks: %d\n", FastVol.NbBlocks);
      for (int block = 1; block <= FastVol.NbBlocks; ++block) {
        printf("NbOfXCells[%d]: %d\n", block, BlkNbXCells[block]);
        printf("NbOfYCells[%d]: %d\n", block, BlkNbYCells[block]);
        printf("NbOfZCells[%d]: %d\n", block, BlkNbZCells[block]);
        printf("LZ[%d]: %le\n", block, BlkLZ[block]);
        printf("CornerZ[%d]: %le\n", block, BlkCrnrZ[block]);
      }
      printf("NbOfOmitVols: %d\n", FastVol.NbOmitVols);
      if (FastVol.NbOmitVols) {
        for (int omit = 1; omit <= FastVol.NbOmitVols; ++omit) {
          printf("OmitVolLX[%d]: %le\n", omit, OmitVolLX[omit]);
          printf("OmitVolLY[%d]: %le\n", omit, OmitVolLY[omit]);
          printf("OmitVolLZ[%d]: %le\n", omit, OmitVolLZ[omit]);
          printf("OmitCrnrX[%d]: %le\n", omit, OmitVolCrnrX[omit]);
          printf("OmitCrnrY[%d]: %le\n", omit, OmitVolCrnrY[omit]);
          printf("OmitCrnrZ[%d]: %le\n", omit, OmitVolCrnrZ[omit]);
        }
      }
      printf("MaxXCells, MaxYCells, MaxZCells: %d, %d, %d\n", MaxXCells,
             MaxYCells, MaxZCells);
    }  // dbgFn
 
    // Following memory allocations are necessary even for creating
    // the fast volumes.
    // In fact, since these tensors are already in place, the fast volumes
    // can be used to evaluate potential and fields immediately after they are
    // created.
    /* Memory wastage!!! Improve as soon as possible. */
    FastPot = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1, MaxYCells + 1,
                       1, MaxZCells + 1);
    FastFX = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1, MaxYCells + 1,
                      1, MaxZCells + 1);
    FastFY = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1, MaxYCells + 1,
                      1, MaxZCells + 1);
    FastFZ = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1, MaxYCells + 1,
                      1, MaxZCells + 1);
 
    if (OptStaggerFastVol) {
      /* Memory wastage!!! Improve as soon as possible. */
      StgFastPot = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1,
                            MaxYCells + 1, 1, MaxZCells + 1);
      StgFastFX = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1,
                           MaxYCells + 1, 1, MaxZCells + 1);
      StgFastFY = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1,
                           MaxYCells + 1, 1, MaxZCells + 1);
      StgFastFZ = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1,
                           MaxYCells + 1, 1, MaxZCells + 1);
    }  // if OptStaggerFastVol
 
    if (OptCreateFastPF) {
      int fstatus = CreateFastVolPF();
      if (fstatus != 0) {
        neBEMMessage("neBEMSolve - Failure computing FastVolPF");
        return -1;
      }
 
      clock_t stopFastClock = clock();
      neBEMTimeElapsed(startFastClock, stopFastClock);
      printf("to compute FastVolPF\n");
    }  // if OptCreateFastPF
 
    if (OptReadFastPF) {  // reading from file
      int nbXCells, nbYCells, nbZCells;
      int tmpblk;
      double xpt, ypt, zpt;
 
      char FastVolPFFile[256];
      strcpy(FastVolPFFile, BCOutDir);
      strcat(FastVolPFFile, "/FastVolPF.out");
      FILE *fFastVolPF = fopen(FastVolPFFile, "r");
      if (fFastVolPF == NULL) {
        neBEMMessage("in neBEMSolve - FastVolPFFile");
        return -1;
      }
 
      fscanf(fFastVolPF, "#block\tX\tY\tZ\tPot\tFX\tFY\tFZ\n");
 
      for (int block = 1; block <= FastVol.NbBlocks; ++block) {
        nbXCells = BlkNbXCells[block];
        nbYCells = BlkNbYCells[block];
        nbZCells = BlkNbZCells[block];
 
        for (int i = 1; i <= nbXCells + 1; ++i) {
          for (int j = 1; j <= nbYCells + 1; ++j) {
            for (int k = 1; k <= nbZCells + 1; ++k) {
              fscanf(fFastVolPF, "%d\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",
                     &tmpblk, &xpt, &ypt, &zpt, &FastPot[block][i][j][k],
                     &FastFX[block][i][j][k], &FastFY[block][i][j][k],
                     &FastFZ[block][i][j][k]);
            }  // loop k
          }    // loop j
        }      // loop i
      }        // loop block
      fclose(fFastVolPF);
 
      if (OptStaggerFastVol) {
        char StgFastVolPFFile[256];
        strcpy(StgFastVolPFFile, BCOutDir);
        strcat(StgFastVolPFFile, "/StgFastVolPF.out");
        FILE *fStgFastVolPF = fopen(StgFastVolPFFile, "r");
        if (fStgFastVolPF == NULL) {
          neBEMMessage("in neBEMSolve - StgFastVolPFFile");
          return -1;
        }
 
        fscanf(fStgFastVolPF, "#block\tX\tY\tZ\tPot\tFX\tFY\tFZ\n");
 
        for (int block = 1; block <= FastVol.NbBlocks; ++block) {
          nbXCells = BlkNbXCells[block];
          nbYCells = BlkNbYCells[block];
          nbZCells = BlkNbZCells[block];
 
          for (int i = 1; i <= nbXCells + 1; ++i) {
            for (int j = 1; j <= nbYCells + 1; ++j) {
              for (int k = 1; k <= nbZCells + 1; ++k) {
                fscanf(fStgFastVolPF, "%d\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",
                       &tmpblk, &xpt, &ypt, &zpt, &StgFastPot[block][i][j][k],
                       &StgFastFX[block][i][j][k], &StgFastFY[block][i][j][k],
                       &StgFastFZ[block][i][j][k]);
              }  // loop k
            }    // loop j
          }      // loop i
        }        // loop block
        fclose(fStgFastVolPF);
      }  // if OptStaggerFastVol
 
      clock_t stopFastClock = clock();
      neBEMTimeElapsed(startFastClock, stopFastClock);
      printf("to read FastVolPF\n");
    }  // if OptReadFastPF
  }    // if OptFastVol
 
  return (0);
}  // neBEMSolve ends

neBEMTimeElapsed()

INTFACEGLOBAL double neBEMTimeElapsed	(	clock_t	t0,
		clock_t	t1 )

Definition at line 3233 of file neBEMInterface.c.

                                                {
  double elapsedTime = ((double)(t1 - t0)) / CLOCKS_PER_SEC;
  printf("neBEMV%s TimeElapsed ===> %lg seconds ", neBEMVersion, elapsedTime);
 
  return (elapsedTime);
}

Referenced by ComputeSolution(), InitialConditions(), InvertMatrix(), neBEMBoundaryInitialConditions(), neBEMDiscretize(), neBEMPrepareWeightingField(), neBEMReadGeometry(), and neBEMSolve().

neBEMVolumeCharge()

INTFACEGLOBAL double neBEMVolumeCharge

(

int

volume

)

Definition at line 2639 of file neBEMInterface.c.

                                     {
  // Initialise the sum
  double sumcharge = 0.0;
 
  // Loop over all elements
  for (int elem = 1; elem <= NbElements; ++elem) {
    // Find the primitive number for the element
    int prim = (EleArr + elem - 1)->PrimitiveNb;
 
    // Find out to which volume this belongs to
    int volref = VolRef1[prim];
 
    // Skip the rest if the volume is not right
    if (volref + 1 != volume) {
      continue;
    }
 
    // Calculate the periodic weight of the primitive
    int rptCnt = 0;
    if (PeriodicInX[prim] || PeriodicInY[prim] || PeriodicInZ[prim]) {
      for (int xrpt = -PeriodicInX[prim]; xrpt <= PeriodicInX[prim]; ++xrpt)
        for (int yrpt = -PeriodicInY[prim]; yrpt <= PeriodicInY[prim]; ++yrpt)
          for (int zrpt = -PeriodicInZ[prim]; zrpt <= PeriodicInZ[prim];
               ++zrpt) {
            if ((xrpt == 0) && (yrpt == 0) && (zrpt == 0))
              continue;
            else
              ++rptCnt;
          }
    } else {
      rptCnt = 1;
    }
 
    // Add the charge
    // printf("Element: %d, volume: %d, charge: %g\n", elem, volref,
    //    (EleArr+elem-1)->Solution * (EleArr+elem-1)->G.dA);
    sumcharge +=
        rptCnt * (EleArr + elem - 1)->Solution * (EleArr + elem - 1)->G.dA;
  }  // loop over elements
 
  // Return the result
  // printf("Charge on volume %d: %g C\n", volume, sumcharge);
  return sumcharge;
}  // end of neBEMVolumeCharge

neBEMVolumeDescription()

INTFACEGLOBAL int neBEMVolumeDescription	(	int	ivol,
		int *	shape,
		int *	material,
		double *	epsilon,
		double *	potential,
		double *	charge,
		int *	boundarytype )

neBEMVolumePoint()

INTFACEGLOBAL int neBEMVolumePoint	(	double	x,
		double	y,
		double	z )

neBEMVolumePrimitives()

INTFACEGLOBAL void neBEMVolumePrimitives	(	int	volume,
		int *	nprim,
		int	primlist[] )

neBEMWeightingField()

INTFACEGLOBAL double neBEMWeightingField	(	Point3D *	point,
		Vector3D *	field,
		int	IdWtField )

Definition at line 2605 of file neBEMInterface.c.

                                                                           {
  double potential;
 
  if (neBEMState < 9) {
    printf("neBEMWeightingField cannot be called before reaching state 9.\n");
    return (-1);
  }
 
  if (OptFixedWtField[IdWtField]) {
    // minimum computation, too restricted! Does not even consider variation
    // through IdWtField
    potential = FixedWtPotential[IdWtField];
    field->X = FixedWtFieldX[IdWtField];
    field->Y = FixedWtFieldY[IdWtField];
    field->Z = FixedWtFieldZ[IdWtField];
  } else if (OptWtFldFastVol[IdWtField]) {
    // bit more computation, lot more flexibility
    // Note: this is not the Create or Read option
    int fstatus = WtFldFastPFAtPoint(point, &potential, field, IdWtField);
    if (fstatus != 0) {
      neBEMMessage("neBEMWeightingField - WtFldFastPFAtPoint");
      return DBL_MAX;
    }
  } else {
    int fstatus = WtFldPFAtPoint(point, &potential, field, IdWtField);
    if (fstatus != 0) {
      neBEMMessage("neBEMWeightingField - WtFldPFAtPoint");
      return DBL_MAX;
    }
  }
 
  return potential;
}  // neBEMWeightingField ends

ReadElements()

INTFACEGLOBAL int ReadElements

(

void

)

Definition at line 3108 of file neBEMInterface.c.

                       {
  char ElementFile[256];
  strcpy(ElementFile, MeshOutDir);
  strcat(ElementFile, "/Elements/StoreElems.out");
 
  FILE *fStrEle = fopen(ElementFile, "r");
  if (fStrEle == NULL) {
    neBEMMessage("ReadElements - Could not open file to read elements");
    return -1;
  }
 
  fscanf(fStrEle, "%d %d\n", &NbSurfs, &NbWires);
 
  for (int prim = 1; prim <= NbPrimitives; ++prim) {
    fscanf(fStrEle, "%d %d\n", &NbSurfSegX[prim], &NbSurfSegZ[prim]);
    fscanf(fStrEle, "%d\n", &NbWireSeg[prim]);
  }
 
  fscanf(fStrEle, "%d\n", &NbElements);
 
  if (neBEMState == 3) {
    if (EleArr) {
      Element *tmp = (Element *)realloc(EleArr, NbElements * sizeof(Element));
      if (tmp != NULL) {
        EleArr = tmp;
        EleCntr = 0;
      } else {
        free(EleArr);
        printf("neBEMDiscretize: Re-allocating EleArr failed.\n");
        fclose(fStrEle);
        return (1);
      }
      printf("neBEMDiscretize: Re-allocated EleArr.\n");
    }  // if EleArr => re-allocation
    else {
      EleArr = (Element *)malloc(NbElements * sizeof(Element));
      if (EleArr == NULL) {
        neBEMMessage("neBEMDiscretize - EleArr malloc");
        return -1;
      }
    }  // else EleArr => fresh allocation
  } else {
    neBEMMessage("neBEMDiscretize - EleArr malloc; neBEMState mismatch!");
    return -1;
  }  // else neBEMState == 3
 
  for (int ele = 1; ele <= NbElements; ++ele) {
    fscanf(fStrEle, "%hd %d %d %d %d\n", &(EleArr + ele - 1)->DeviceNb,
           &(EleArr + ele - 1)->ComponentNb, &(EleArr + ele - 1)->PrimitiveNb,
           &(EleArr + ele - 1)->InterfaceId, &(EleArr + ele - 1)->Id);
    fscanf(fStrEle, "%hd %le %le %le %le %le %le\n",
           &(EleArr + ele - 1)->G.Type, &(EleArr + ele - 1)->G.Origin.X,
           &(EleArr + ele - 1)->G.Origin.Y, &(EleArr + ele - 1)->G.Origin.Z,
           &(EleArr + ele - 1)->G.LX, &(EleArr + ele - 1)->G.LZ,
           &(EleArr + ele - 1)->G.dA);
    fscanf(fStrEle, "%le %le %le\n", &(EleArr + ele - 1)->G.DC.XUnit.X,
           &(EleArr + ele - 1)->G.DC.XUnit.Y,
           &(EleArr + ele - 1)->G.DC.XUnit.Z);
    fscanf(fStrEle, "%le %le %le\n", &(EleArr + ele - 1)->G.DC.YUnit.X,
           &(EleArr + ele - 1)->G.DC.YUnit.Y,
           &(EleArr + ele - 1)->G.DC.YUnit.Z);
    fscanf(fStrEle, "%le %le %le\n", &(EleArr + ele - 1)->G.DC.ZUnit.X,
           &(EleArr + ele - 1)->G.DC.ZUnit.Y,
           &(EleArr + ele - 1)->G.DC.ZUnit.Z);
    fscanf(fStrEle, "%hd %le\n", &(EleArr + ele - 1)->E.Type,
           &(EleArr + ele - 1)->E.Lambda);
    fscanf(fStrEle, "%hd %le %le %le %le\n", &(EleArr + ele - 1)->BC.NbOfBCs,
           &(EleArr + ele - 1)->BC.CollPt.X, &(EleArr + ele - 1)->BC.CollPt.Y,
           &(EleArr + ele - 1)->BC.CollPt.Z, &(EleArr + ele - 1)->BC.Value);
    fscanf(fStrEle, "%le %le\n", &(EleArr + ele - 1)->Solution,
           &(EleArr + ele - 1)->Assigned);
  }
 
  fscanf(fStrEle, "%d %d %d %d\n", &NbPointsKnCh, &NbLinesKnCh, &NbAreasKnCh,
         &NbVolumesKnCh);
 
  for (int pt = 1; pt <= NbPointsKnCh; ++pt) {
    fscanf(fStrEle, "%d %le\n", &(PointKnChArr + pt - 1)->Nb,
           &(PointKnChArr + pt - 1)->Assigned);
    fscanf(fStrEle, "%le %le %le\n", &(PointKnChArr + pt - 1)->P.X,
           &(PointKnChArr + pt - 1)->P.Y, &(PointKnChArr + pt - 1)->P.Z);
  }
 
  for (int line = 1; line <= NbLinesKnCh; ++line) {
    fscanf(fStrEle, "%d %le %le\n", &(LineKnChArr + line - 1)->Nb,
           &(LineKnChArr + line - 1)->Radius,
           &(LineKnChArr + line - 1)->Assigned);
    fscanf(fStrEle, "%le %le %le\n", &(LineKnChArr + line - 1)->Start.X,
           &(LineKnChArr + line - 1)->Start.Y,
           &(LineKnChArr + line - 1)->Start.Z);
    fscanf(fStrEle, "%le %le %le\n", &(LineKnChArr + line - 1)->Stop.X,
           &(LineKnChArr + line - 1)->Stop.Y,
           &(LineKnChArr + line - 1)->Stop.Z);
  }
 
  for (int area = 1; area <= NbAreasKnCh; ++area) {
    fscanf(fStrEle, "%d %d %le\n", &(AreaKnChArr + area - 1)->Nb,
           &(AreaKnChArr + area - 1)->NbVertices,
           &(AreaKnChArr + area - 1)->Assigned);
    for (int vert = 1; vert <= (AreaKnChArr + area - 1)->NbVertices; ++vert) {
      fscanf(fStrEle, "%le %le %le\n",
             &(AreaKnChArr + area - 1)->Vertex[vert].X,
             &(AreaKnChArr + area - 1)->Vertex[vert].Y,
             &(AreaKnChArr + area - 1)->Vertex[vert].Z);
    }
  }
 
  for (int vol = 1; vol <= NbVolumesKnCh; ++vol) {
    fscanf(fStrEle, "%d %d %le\n", &(VolumeKnChArr + vol - 1)->Nb,
           &(VolumeKnChArr + vol - 1)->NbVertices,
           &(VolumeKnChArr + vol - 1)->Assigned);
    for (int vert = 1; vert <= (VolumeKnChArr + vol - 1)->NbVertices; ++vert) {
      fscanf(fStrEle, "%le %le %le\n",
             &(VolumeKnChArr + vol - 1)->Vertex[vert].X,
             &(VolumeKnChArr + vol - 1)->Vertex[vert].Y,
             &(VolumeKnChArr + vol - 1)->Vertex[vert].Z);
    }
  }
 
  fclose(fStrEle);
 
  return 0;
}  // ReadElements ends

Referenced by neBEMDiscretize().

ReadInitFile()

INTFACEGLOBAL int ReadInitFile

(

char

filename[]

)

ReadPrimitives()

INTFACEGLOBAL int ReadPrimitives

(

void

)

Definition at line 3004 of file neBEMInterface.c.

                         {
  int dbgFn = 0;
 
  char PrimitiveFile[256];
 
  strcpy(PrimitiveFile, ModelOutDir);
  strcat(PrimitiveFile, "/Primitives/StorePrims.out");
 
  FILE *fStrPrm = fopen(PrimitiveFile, "r");
  if (fStrPrm == NULL) {
    neBEMMessage("ReadPrimitives - Could not open file to read primitives");
    return -1;
  }
 
  fscanf(fStrPrm, "%d %d\n", &NbVolumes, &VolMax);
  fscanf(fStrPrm, "%d\n", &NbPrimitives);
  fscanf(fStrPrm, "%d\n", &MaxNbVertices);
 
  // assign neBEMState and allocate memory
  neBEMState = 2;
  PrimType = ivector(1, NbPrimitives);
  NbVertices = ivector(1, NbPrimitives);
  XVertex = dmatrix(1, NbPrimitives, 0, MaxNbVertices - 1);
  YVertex = dmatrix(1, NbPrimitives, 0, MaxNbVertices - 1);
  ZVertex = dmatrix(1, NbPrimitives, 0, MaxNbVertices - 1);
  XNorm = dvector(1, NbPrimitives);
  YNorm = dvector(1, NbPrimitives);
  ZNorm = dvector(1, NbPrimitives);
  Radius = dvector(1, NbPrimitives);  // can lead to a little memory misuse
  VolRef1 = ivector(1, NbPrimitives);
  VolRef2 = ivector(1, NbPrimitives);
  NbSurfSegX = ivector(1, NbPrimitives);
  NbSurfSegZ = ivector(1, NbPrimitives);
  NbWireSeg = ivector(1, NbPrimitives);  // little memory misuse
  InterfaceType = ivector(1, NbPrimitives);
  Lambda = dvector(1, NbPrimitives);
  ApplPot = dvector(1, NbPrimitives);
  ApplCh = dvector(1, NbPrimitives);
  PeriodicTypeX = ivector(1, NbPrimitives);
  PeriodicTypeY = ivector(1, NbPrimitives);
  PeriodicTypeZ = ivector(1, NbPrimitives);
  PeriodicInX = ivector(1, NbPrimitives);
  PeriodicInY = ivector(1, NbPrimitives);
  PeriodicInZ = ivector(1, NbPrimitives);
  XPeriod = dvector(1, NbPrimitives);
  YPeriod = dvector(1, NbPrimitives);
  ZPeriod = dvector(1, NbPrimitives);
  MirrorDistXFromOrigin = dvector(1, NbPrimitives);
  MirrorDistYFromOrigin = dvector(1, NbPrimitives);
  MirrorDistZFromOrigin = dvector(1, NbPrimitives);
 
  for (int prim = 1; prim <= NbPrimitives; ++prim) {
    fscanf(fStrPrm, "%d\n", &PrimType[prim]);
    fscanf(fStrPrm, "%d\n", &InterfaceType[prim]);
    fscanf(fStrPrm, "%d\n", &NbVertices[prim]);
 
    for (int vert = 0; vert < NbVertices[prim]; ++vert) {
      fscanf(fStrPrm, "%le %le %le\n", &XVertex[prim][vert],
             &YVertex[prim][vert], &ZVertex[prim][vert]);
    }  // vert loop
 
    fscanf(fStrPrm, "%le %le %le\n", &XNorm[prim], &YNorm[prim], &ZNorm[prim]);
    fscanf(fStrPrm, "%le\n", &Radius[prim]);
 
    fscanf(fStrPrm, "%le %le %le %le %le\n", &Epsilon1[prim], &Epsilon2[prim],
           &Lambda[prim], &ApplPot[prim], &ApplCh[prim]);
 
    fscanf(fStrPrm, "%d %d\n", &VolRef1[prim], &VolRef2[prim]);
 
    fscanf(fStrPrm, "%d %d %d\n", &PeriodicTypeX[prim], &PeriodicTypeY[prim],
           &PeriodicTypeZ[prim]);
    fscanf(fStrPrm, "%d %d %d\n", &PeriodicInX[prim], &PeriodicInY[prim],
           &PeriodicInZ[prim]);
    fscanf(fStrPrm, "%le %le %le\n", &XPeriod[prim], &YPeriod[prim],
           &ZPeriod[prim]);
    fscanf(fStrPrm, "%le %le %le\n", &MirrorDistXFromOrigin[prim],
           &MirrorDistYFromOrigin[prim], &MirrorDistZFromOrigin[prim]);
  }  // prim loop
 
  volRef = ivector(0, VolMax);
  volShape = ivector(0, VolMax);
  volMaterial = ivector(0, VolMax);
  volEpsilon = dvector(0, VolMax);
  volPotential = dvector(0, VolMax);
  volCharge = dvector(0, VolMax);
  volBoundaryType = ivector(0, VolMax);
  for (int volref = 0; volref <= VolMax; ++volref) {
    neBEMVolumeDescription(volref, &volShape[volref], &volMaterial[volref],
                           &volEpsilon[volref], &volPotential[volref],
                           &volCharge[volref], &volBoundaryType[volref]);
    if (dbgFn) {
      printf("volref: %d\n", volref);
      printf("shape: %d,  material: %d\n", volShape[volref],
             volMaterial[volref]);
      printf("eps: %lg,  pot: %lg\n", volEpsilon[volref], volPotential[volref]);
      printf("q: %lg,  type: %d\n", volCharge[volref], volBoundaryType[volref]);
    }
  }  // volume loop
 
  fclose(fStrPrm);
 
  return 0;
}  // ReadPrimitives ends

Referenced by neBEMReadGeometry().

WriteElements()

INTFACEGLOBAL int WriteElements

(

void

)

Definition at line 2910 of file neBEMInterface.c.

                        {
  char ElementFile[256];
 
  strcpy(ElementFile, MeshOutDir);
  strcat(ElementFile, "/Elements/StoreElems.out");
 
  FILE *fStrEle = fopen(ElementFile, "w");
  if (fStrEle == NULL) {
    neBEMMessage("WriteElements - Could not create file to store elements");
    return -1;
  }
 
  fprintf(fStrEle, "%d %d\n", NbSurfs, NbWires);
 
  for (int prim = 1; prim <= NbPrimitives; ++prim) {
    fprintf(fStrEle, "%d %d\n", NbSurfSegX[prim], NbSurfSegZ[prim]);
    fprintf(fStrEle, "%d\n", NbWireSeg[prim]);
  }
 
  fprintf(fStrEle, "%d\n", NbElements);
 
  for (int ele = 1; ele <= NbElements; ++ele) {
    fprintf(fStrEle, "%d %d %d %d %d\n", (EleArr + ele - 1)->DeviceNb,
            (EleArr + ele - 1)->ComponentNb, (EleArr + ele - 1)->PrimitiveNb,
            (EleArr + ele - 1)->InterfaceId, (EleArr + ele - 1)->Id);
    fprintf(fStrEle, "%d %le %le %le %le %le %le\n", (EleArr + ele - 1)->G.Type,
            (EleArr + ele - 1)->G.Origin.X, (EleArr + ele - 1)->G.Origin.Y,
            (EleArr + ele - 1)->G.Origin.Z, (EleArr + ele - 1)->G.LX,
            (EleArr + ele - 1)->G.LZ, (EleArr + ele - 1)->G.dA);
    fprintf(fStrEle, "%le %le %le\n", (EleArr + ele - 1)->G.DC.XUnit.X,
            (EleArr + ele - 1)->G.DC.XUnit.Y, (EleArr + ele - 1)->G.DC.XUnit.Z);
    fprintf(fStrEle, "%le %le %le\n", (EleArr + ele - 1)->G.DC.YUnit.X,
            (EleArr + ele - 1)->G.DC.YUnit.Y, (EleArr + ele - 1)->G.DC.YUnit.Z);
    fprintf(fStrEle, "%le %le %le\n", (EleArr + ele - 1)->G.DC.ZUnit.X,
            (EleArr + ele - 1)->G.DC.ZUnit.Y, (EleArr + ele - 1)->G.DC.ZUnit.Z);
    fprintf(fStrEle, "%d %le\n", (EleArr + ele - 1)->E.Type,
            (EleArr + ele - 1)->E.Lambda);
    fprintf(fStrEle, "%d %le %le %le %le\n", (EleArr + ele - 1)->BC.NbOfBCs,
            (EleArr + ele - 1)->BC.CollPt.X, (EleArr + ele - 1)->BC.CollPt.Y,
            (EleArr + ele - 1)->BC.CollPt.Z, (EleArr + ele - 1)->BC.Value);
    fprintf(fStrEle, "%le %le\n", (EleArr + ele - 1)->Solution,
            (EleArr + ele - 1)->Assigned);
  }
 
  fprintf(fStrEle, "%d %d %d %d\n", NbPointsKnCh, NbLinesKnCh, NbAreasKnCh,
          NbVolumesKnCh);
 
  for (int pt = 1; pt <= NbPointsKnCh; ++pt) {
    fprintf(fStrEle, "%d %le\n", (PointKnChArr + pt - 1)->Nb,
            (PointKnChArr + pt - 1)->Assigned);
    fprintf(fStrEle, "%le %le %le\n", (PointKnChArr + pt - 1)->P.X,
            (PointKnChArr + pt - 1)->P.Y, (PointKnChArr + pt - 1)->P.Z);
  }
 
  for (int line = 1; line <= NbLinesKnCh; ++line) {
    fprintf(fStrEle, "%d %le %le\n", (LineKnChArr + line - 1)->Nb,
            (LineKnChArr + line - 1)->Radius,
            (LineKnChArr + line - 1)->Assigned);
    fprintf(fStrEle, "%le %le %le\n", (LineKnChArr + line - 1)->Start.X,
            (LineKnChArr + line - 1)->Start.Y,
            (LineKnChArr + line - 1)->Start.Z);
    fprintf(fStrEle, "%le %le %le\n", (LineKnChArr + line - 1)->Stop.X,
            (LineKnChArr + line - 1)->Stop.Y, (LineKnChArr + line - 1)->Stop.Z);
  }
 
  for (int area = 1; area <= NbAreasKnCh; ++area) {
    fprintf(fStrEle, "%d %d %le\n", (AreaKnChArr + area - 1)->Nb,
            (AreaKnChArr + area - 1)->NbVertices,
            (AreaKnChArr + area - 1)->Assigned);
    for (int vert = 1; vert <= (AreaKnChArr + area - 1)->NbVertices; ++vert) {
      fprintf(fStrEle, "%le %le %le\n",
              (AreaKnChArr + area - 1)->Vertex[vert].X,
              (AreaKnChArr + area - 1)->Vertex[vert].Y,
              (AreaKnChArr + area - 1)->Vertex[vert].Z);
    }
  }
 
  for (int vol = 1; vol <= NbVolumesKnCh; ++vol) {
    fprintf(fStrEle, "%d %d %le\n", (VolumeKnChArr + vol - 1)->Nb,
            (VolumeKnChArr + vol - 1)->NbVertices,
            (VolumeKnChArr + vol - 1)->Assigned);
    for (int vert = 1; vert <= (VolumeKnChArr + vol - 1)->NbVertices; ++vert) {
      fprintf(fStrEle, "%le %le %le\n",
              (VolumeKnChArr + vol - 1)->Vertex[vert].X,
              (VolumeKnChArr + vol - 1)->Vertex[vert].Y,
              (VolumeKnChArr + vol - 1)->Vertex[vert].Z);
    }
  }
 
  fclose(fStrEle);
 
  return 0;
}  // WriteElements ends

Referenced by neBEMDiscretize().

WritePrimitives()

INTFACEGLOBAL int WritePrimitives

(

void

)

Definition at line 2861 of file neBEMInterface.c.

                          {
  char PrimitiveFile[256];
 
  strcpy(PrimitiveFile, ModelOutDir);
  strcat(PrimitiveFile, "/Primitives/StorePrims.out");
 
  FILE *fStrPrm = fopen(PrimitiveFile, "w");
  if (fStrPrm == NULL) {
    neBEMMessage("WritePrimitives - Could not create file to store primitives");
    return -1;
  }
 
  fprintf(fStrPrm, "%d %d\n", NbVolumes, VolMax);
  fprintf(fStrPrm, "%d\n", NbPrimitives);
  fprintf(fStrPrm, "%d\n", MaxNbVertices);
 
  for (int prim = 1; prim <= NbPrimitives; ++prim) {
    fprintf(fStrPrm, "%d\n", PrimType[prim]);
    fprintf(fStrPrm, "%d\n", InterfaceType[prim]);
    fprintf(fStrPrm, "%d\n", NbVertices[prim]);
 
    for (int vert = 0; vert < NbVertices[prim]; ++vert) {
      fprintf(fStrPrm, "%le %le %le\n", XVertex[prim][vert],
              YVertex[prim][vert], ZVertex[prim][vert]);
    }  // vert loop
 
    fprintf(fStrPrm, "%le %le %le\n", XNorm[prim], YNorm[prim], ZNorm[prim]);
    fprintf(fStrPrm, "%le\n", Radius[prim]);
 
    fprintf(fStrPrm, "%le %le %le %le %le\n", Epsilon1[prim], Epsilon2[prim],
            Lambda[prim], ApplPot[prim], ApplCh[prim]);
 
    fprintf(fStrPrm, "%d %d\n", VolRef1[prim], VolRef2[prim]);
 
    fprintf(fStrPrm, "%d %d %d\n", PeriodicTypeX[prim], PeriodicTypeY[prim],
            PeriodicTypeZ[prim]);
    fprintf(fStrPrm, "%d %d %d\n", PeriodicInX[prim], PeriodicInY[prim],
            PeriodicInZ[prim]);
    fprintf(fStrPrm, "%le %le %le\n", XPeriod[prim], YPeriod[prim],
            ZPeriod[prim]);
    fprintf(fStrPrm, "%le %le %le\n", MirrorDistXFromOrigin[prim],
            MirrorDistYFromOrigin[prim], MirrorDistZFromOrigin[prim]);
  }  // prim loop
 
  fclose(fStrPrm);
 
  return 0;
}  // WritePrimitives ends

Referenced by neBEMReadGeometry().

Variable Documentation

DeviceInputFile

INTFACEGLOBAL char DeviceInputFile[256]

Definition at line 45 of file neBEMInterface.h.

Referenced by neBEMInitialize(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

fgnuElem

INTFACEGLOBAL FILE * fgnuElem

Definition at line 55 of file neBEMInterface.h.

Referenced by DiscretizeRectangle(), DiscretizeTriangle(), and neBEMDiscretize().

fgnuMesh

INTFACEGLOBAL FILE * fgnuMesh

Definition at line 55 of file neBEMInterface.h.

Referenced by DiscretizeRectangle(), DiscretizeTriangle(), and neBEMDiscretize().

fgnuPrim

INTFACEGLOBAL FILE* fgnuPrim

Definition at line 55 of file neBEMInterface.h.

Referenced by DiscretizeRectangle(), DiscretizeTriangle(), DiscretizeWire(), and neBEMDiscretize().

NbThreads

INTFACEGLOBAL int NbThreads

Definition at line 49 of file neBEMInterface.h.

Referenced by neBEMInitialize().

neBEMPFCallCntr

INTFACEGLOBAL int neBEMPFCallCntr

Definition at line 109 of file neBEMInterface.h.

neBEMState

INTFACEGLOBAL int neBEMState

Definition at line 23 of file neBEMInterface.h.

Referenced by DecomposeMatrixSVD(), InvertMatrix(), LHMatrix(), neBEMBoundaryInitialConditions(), neBEMDiscretize(), neBEMInitialize(), neBEMPF(), neBEMPrepareWeightingField(), neBEMReadGeometry(), neBEM::neBEMSetDefaults(), neBEMSolve(), neBEMWeightingField(), ReadElements(), ReadInvertedMatrix(), ReadPrimitives(), ReadSolution(), RHVector(), and Solve().

OptDeviceFile

INTFACEGLOBAL int OptDeviceFile

Definition at line 44 of file neBEMInterface.h.

Referenced by neBEMInitialize(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

OptElementFiles

INTFACEGLOBAL int OptElementFiles

Definition at line 98 of file neBEMInterface.h.

Referenced by DiscretizeRectangle(), DiscretizeTriangle(), DiscretizeWire(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

OptGnuplot

INTFACEGLOBAL int OptGnuplot

Definition at line 54 of file neBEMInterface.h.

Referenced by DiscretizeRectangle(), DiscretizeTriangle(), DiscretizeWire(), neBEMDiscretize(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

OptGnuplotElements

INTFACEGLOBAL int OptGnuplotElements

Definition at line 54 of file neBEMInterface.h.

Referenced by DiscretizeRectangle(), DiscretizeTriangle(), DiscretizeWire(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

OptGnuplotPrimitives

INTFACEGLOBAL int OptGnuplotPrimitives

Definition at line 54 of file neBEMInterface.h.

Referenced by DiscretizeRectangle(), DiscretizeTriangle(), DiscretizeWire(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

OptPrimitiveFiles

INTFACEGLOBAL int OptPrimitiveFiles

Definition at line 62 of file neBEMInterface.h.

Referenced by DiscretizeRectangle(), DiscretizeTriangle(), DiscretizeWire(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

OptPrintPrimaryDetails

INTFACEGLOBAL int OptPrintPrimaryDetails

Definition at line 52 of file neBEMInterface.h.

Referenced by neBEMReadGeometry(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

OptPrintVertexAndNormal

INTFACEGLOBAL int OptPrintVertexAndNormal

Definition at line 62 of file neBEMInterface.h.

Referenced by DiscretizeRectangle(), DiscretizeTriangle(), DiscretizeWire(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

OptPrintVolumeDetails

INTFACEGLOBAL int OptPrintVolumeDetails

Definition at line 52 of file neBEMInterface.h.

Referenced by neBEMReadGeometry(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

OptReuseDir

INTFACEGLOBAL int OptReuseDir

Definition at line 142 of file neBEMInterface.h.

Referenced by CreateDirStr(), neBEM::neBEMSetDefaults(), and neBEM::ReadInitFile().

startClock

INTFACEGLOBAL clock_t startClock

Definition at line 19 of file neBEMInterface.h.

Referenced by ComputeSolution(), InitialConditions(), neBEMBoundaryInitialConditions(), neBEMDiscretize(), and neBEMReadGeometry().

stopClock

INTFACEGLOBAL clock_t stopClock

Definition at line 19 of file neBEMInterface.h.

Referenced by ComputeSolution(), InitialConditions(), neBEMBoundaryInitialConditions(), neBEMDiscretize(), and neBEMReadGeometry().