Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
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)   (sqrt(p.X * p.X + p.Y * p.Y + p.Z * p.Z))
 

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 neBEMBoundaryConditions (void)
 
INTFACEGLOBAL int neBEMKnownCharges (void)
 
INTFACEGLOBAL int neBEMChargingUp (int InfluenceMatrixFlag)
 
INTFACEGLOBAL int neBEMSolve (void)
 
INTFACEGLOBAL int neBEMField (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 neBEMFieldCallCntr
 
INTFACEGLOBAL int OptReuseDir
 

Macro Definition Documentation

◆ INTFACEGLOBAL

#define INTFACEGLOBAL   extern

Definition at line 7 of file neBEMInterface.h.

◆ MODULUS

#define MODULUS (   p)    (sqrt(p.X * p.X + p.Y * p.Y + p.Z * p.Z))

Definition at line 35 of file neBEMInterface.h.

◆ 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.

Function Documentation

◆ CreateDirStr()

INTFACEGLOBAL int CreateDirStr ( void  )

Definition at line 4774 of file neBEMInterface.c.

4774 {
4775 char strModelCntr[10], strMeshCntr[10], strBCCntr[10], strPPCntr[10];
4776 int CreateOrUseDir(char[]);
4777 int CreateDirOrQuit(char[]);
4778
4779 sprintf(strModelCntr, "/Model%d", ModelCntr);
4780 sprintf(strMeshCntr, "/M%d", MeshCntr);
4781 sprintf(strBCCntr, "/BC%d", BCCntr);
4782 sprintf(strPPCntr, "/PP%d", PPCntr);
4783
4784 strcpy(ModelOutDir, DeviceOutDir);
4785 strcat(ModelOutDir, strModelCntr);
4786 strcpy(NativeOutDir, ModelOutDir);
4787 strcat(NativeOutDir, "/neBEMNatives/");
4788 strcpy(NativePrimDir, NativeOutDir);
4789 strcat(NativePrimDir, "Primitives/");
4790 strcpy(MeshOutDir, ModelOutDir);
4791 strcat(MeshOutDir, strMeshCntr);
4792 strcpy(BCOutDir, MeshOutDir);
4793 strcat(BCOutDir, strBCCntr);
4794 strcpy(PPOutDir, BCOutDir);
4795 strcat(PPOutDir, strPPCntr);
4796
4797 // Create DeviceOutDir, if necessary
4798 int fstatus = CreateOrUseDir(DeviceOutDir);
4799 if (fstatus != 0) {
4800 neBEMMessage("CreateDirStr - CreateOrUseDir");
4801 return -1;
4802 }
4803
4804 if (NewModel) {
4805 // create ModelOutDir
4806 if (OptReuseDir) {
4807 fstatus = CreateOrUseDir(ModelOutDir);
4808 fstatus = CreateOrUseDir(NativeOutDir);
4809 fstatus = CreateOrUseDir(NativePrimDir);
4810 } else {
4811 fstatus = CreateDirOrQuit(ModelOutDir);
4812 fstatus = CreateDirOrQuit(NativeOutDir);
4813 fstatus = CreateDirOrQuit(NativePrimDir);
4814 }
4815 if (fstatus != 0) {
4816 neBEMMessage("CreateDirStr - ModelOutDir");
4817 return -1;
4818 }
4819 }
4820
4821 if (NewMesh) {
4822 // create MeshOutDir
4823 if (OptReuseDir)
4824 fstatus = CreateOrUseDir(MeshOutDir);
4825 else
4826 fstatus = CreateDirOrQuit(MeshOutDir);
4827 if (fstatus != 0) {
4828 neBEMMessage("CreateDirStr - MeshOutDir");
4829 return -1;
4830 }
4831 }
4832
4833 if (NewBC) {
4834 // create BCOutDir
4835 if (OptReuseDir)
4836 fstatus = CreateOrUseDir(BCOutDir);
4837 else
4838 fstatus = CreateDirOrQuit(BCOutDir);
4839 if (fstatus != 0) {
4840 neBEMMessage("CreateDirStr - BCOutDir");
4841 return -1;
4842 }
4843 }
4844
4845 if (NewPP) {
4846 // create PPOutDir
4847 if (OptReuseDir)
4848 fstatus = CreateOrUseDir(PPOutDir);
4849 else
4850 fstatus = CreateDirOrQuit(PPOutDir);
4851 if (fstatus != 0) {
4852 neBEMMessage("CreateDirStr - PPOutDir");
4853 return -1;
4854 }
4855 }
4856
4857 // Create other relevant sub-directories
4858 char subdir[256];
4859
4860 strcpy(subdir, ModelOutDir);
4861 strcat(subdir, "/Primitives/");
4862 if (OptReuseDir)
4863 fstatus = CreateOrUseDir(subdir);
4864 else
4865 fstatus = CreateDirOrQuit(subdir);
4866
4867 strcpy(subdir, MeshOutDir);
4868 strcat(subdir, "/Elements/");
4869 if (OptReuseDir)
4870 fstatus = CreateOrUseDir(subdir);
4871 else
4872 fstatus = CreateDirOrQuit(subdir);
4873
4874 strcpy(subdir, MeshOutDir);
4875 strcat(subdir, "/GViewDir/");
4876 if (OptReuseDir)
4877 fstatus = CreateOrUseDir(subdir);
4878 else
4879 fstatus = CreateDirOrQuit(subdir);
4880
4881 return (0);
4882} // CreateDirStr ends
int neBEMMessage(const char *message)
int CreateOrUseDir(char dirname[])
int CreateDirOrQuit(char dirname[])
INTFACEGLOBAL int OptReuseDir
neBEMGLOBAL int MeshCntr
Definition: neBEM.h:224
neBEMGLOBAL int PPCntr
Definition: neBEM.h:224
neBEMGLOBAL char MeshOutDir[256]
Definition: neBEM.h:263
neBEMGLOBAL int NewModel
Definition: neBEM.h:220
neBEMGLOBAL char PPOutDir[256]
Definition: neBEM.h:263
neBEMGLOBAL int NewMesh
Definition: neBEM.h:220
neBEMGLOBAL char ModelOutDir[256]
Definition: neBEM.h:262
neBEMGLOBAL int NewPP
Definition: neBEM.h:220
neBEMGLOBAL char NativeOutDir[256]
Definition: neBEM.h:262
neBEMGLOBAL int BCCntr
Definition: neBEM.h:224
neBEMGLOBAL int ModelCntr
Definition: neBEM.h:224
neBEMGLOBAL char NativePrimDir[256]
Definition: neBEM.h:263
neBEMGLOBAL char DeviceOutDir[256]
Definition: neBEM.h:262
neBEMGLOBAL char BCOutDir[256]
Definition: neBEM.h:263
neBEMGLOBAL int NewBC
Definition: neBEM.h:220

Referenced by neBEMInitialize().

◆ neBEMBoundaryConditions()

INTFACEGLOBAL int neBEMBoundaryConditions ( void  )

Definition at line 2214 of file neBEMInterface.c.

2214 {
2215 startClock = clock();
2216
2217 // The boundary conditions can be set only with neBEMState == 4 or 7
2218 // This state is assigned either after element discretization has been
2219 // completed or in a condition when we are looking for modifying only the
2220 // boundary condition for a device having same geometry (hence, the same
2221 // inverted influence coefficient matrix)
2222 if ((neBEMState == 4) || (neBEMState == 7)) {
2223 int fstatus = BoundaryConditions();
2224 if (fstatus != 0) {
2225 neBEMMessage("neBEMBondaryConditions - BoundaryConditions");
2226 return -1;
2227 }
2228 if (neBEMState == 4) neBEMState = 5; // create LHMatrix, invert etc
2229 if (neBEMState == 7) neBEMState = 8; // assume LHS and inversion to be over
2230 } else {
2231 printf("Boundary conditions can be set only in state 4 / 7 ...\n");
2232 printf("returning ...\n");
2233 return (-1);
2234 }
2235
2236 stopClock = clock();
2238 printf("to setup boundary conditions.\n");
2239
2240 return (0);
2241} // neBEMBoundaryConditions ends
double neBEMTimeElapsed(clock_t t0, clock_t t1)
INTFACEGLOBAL clock_t stopClock
INTFACEGLOBAL int neBEMState
INTFACEGLOBAL clock_t startClock
neBEMGLOBAL int BoundaryConditions(void)
Definition: ReTriM.c:2121

◆ neBEMChargingUp()

INTFACEGLOBAL int neBEMChargingUp ( int  InfluenceMatrixFlag)

Definition at line 2531 of file neBEMInterface.c.

2531 {
2532 int debugFn = 0;
2533
2534 // status of elements before being charged up
2535 if (debugFn) {
2536 for (int ele = 1; ele <= NbElements; ++ele) {
2537 printf("ele, Assigned charge: %d, %lg\n", ele,
2538 (EleArr + ele - 1)->Assigned);
2539 }
2540 }
2541
2542 // Set up parameters related to charging-up calculations
2543 // The plan is to distribute the electrons and ions ending in dielectric
2544 // volumes to the elements on the volumes
2545 // Only the electrons, to begin with
2546 {
2547 FILE *ChargingUpInpFile = fopen("neBEMChargingUp.inp", "r");
2548 if (ChargingUpInpFile == NULL) {
2549 printf(
2550 "neBEMChargingUp.inp absent ... assuming no charging up effect "
2551 "...\n");
2552 // assign NbChUpEEle and NbChUpIEle and Prims to be zeros?
2553 } else {
2554 fscanf(ChargingUpInpFile, "OptChargingUp: %d\n", &OptChargingUp);
2555 if (!OptChargingUp)
2556 printf("OptChargingUp = 0 ... assuming no charging up effect ...\n");
2557 if (1) printf("OptChargingUp: %d\n", OptChargingUp);
2558
2559 if (OptChargingUp) {
2560 char ChargingUpEFile[256];
2561 char ChargingUpIFile[256];
2562 double ChUpFactor;
2563 int *NbChUpEonEle, *NbChUpIonEle;
2564
2565 fscanf(ChargingUpInpFile, "ChargingUpEFile: %255s\n", ChargingUpEFile);
2566 fscanf(ChargingUpInpFile, "ChargingUpIFile: %255s\n", ChargingUpIFile);
2567 fscanf(ChargingUpInpFile, "ChUpFactor: %lg\n", &ChUpFactor);
2568 if (1) {
2569 printf("ChargingUpEFile: %s\n", ChargingUpEFile);
2570 printf("ChargingUpIFile: %s\n", ChargingUpIFile);
2571 printf("ChUpFactor: %lg\n", ChUpFactor);
2572 }
2573
2574 { // Calculation for electrons
2575 FILE *fptrChargingUpEFile = fopen(ChargingUpEFile, "r");
2576 if (fptrChargingUpEFile == NULL) {
2577 neBEMMessage("ChargingUpE file absent ... returning\n");
2578 return -10;
2579 }
2580 int NbOfE = neBEMGetNbOfLines(ChargingUpEFile);
2581 if (NbOfE <= 1) {
2582 neBEMMessage("Too few lines in ChargingUpE ... returning\n");
2583 fclose(fptrChargingUpEFile);
2584 return -11;
2585 }
2586 NbChUpEonEle = (int *)malloc((NbElements + 1) * sizeof(int));
2587 if (NbChUpEonEle == NULL) {
2588 neBEMMessage("Memory allocation failed ... returning\n");
2589 fclose(fptrChargingUpEFile);
2590 return -11;
2591 }
2592 for (int ele = 0; ele <= NbElements; ++ele) {
2593 // CHECK!!! ele limits start from 0, but all else from 1 to ...
2594 NbChUpEonEle[ele] = 0;
2595 }
2596
2597 // read the header line
2598 char header[256];
2599 fgets(header, 256, fptrChargingUpEFile);
2600
2601 --NbOfE; // one line was for the header
2602 if (debugFn) printf("No. of electrons: %d\n", NbOfE);
2603 const char *tmpEFile = "/tmp/ElectronTempFile.out";
2604 FILE *ftmpEF = fopen(tmpEFile, "w");
2605 if (ftmpEF == NULL) {
2606 printf("cannot open temporary output file ... returning ...\n");
2607 free(NbChUpEonEle);
2608 return -100;
2609 }
2610 FILE *fPtEChUpMap = fopen("PtEChUpMap.out", "w");
2611 if (fPtEChUpMap == NULL) {
2612 printf("cannot open PtEChUpMap.out file for writing ...\n");
2613 free(NbChUpEonEle);
2614 fclose(ftmpEF);
2615 return 110;
2616 }
2617
2618 char label;
2619 int vol, enb; // label, volume and electron number
2620 double xlbend, ylbend, zlbend; // lbend == Last But END
2621 double xend, yend, zend;
2622 Point3D
2623 ptintsct; // each electron likely to have an intersection point
2624 for (int electron = 1; electron <= NbOfE; ++electron) {
2625 fscanf(fptrChargingUpEFile, "%c %d %d %lg %lg %lg %lg %lg %lg\n",
2626 &label, &vol, &enb, &xlbend, &xend, &ylbend, &yend, &zlbend,
2627 &zend);
2628 xlbend *= 0.01;
2629 xend *= 0.01;
2630 ylbend *= 0.01;
2631 yend *= 0.01;
2632 zlbend *= 0.01;
2633 zend *= 0.01;
2634 ptintsct.X = 0.0;
2635 ptintsct.Y = 0.0;
2636 ptintsct.Z = 0.0; // initialize
2637
2638 // find the parametric equation of this last segment
2639 // if xend < xlbend, swap the directions
2640 // This has not been mentioned as mandatory in Vince's book
2641 // "Geometry for Computer Graphics", but implied in the book "A
2642 // Programmer's Geometry"
2643 if (xend < xlbend) // not implemented now
2644 {
2645 }
2646 double lseg = (xend - xlbend) * (xend - xlbend) +
2647 (yend - ylbend) * (yend - ylbend) +
2648 (zend - zlbend) * (zend - zlbend);
2649 lseg = sqrt(lseg);
2650 double xgrd =
2651 (xend - xlbend) / lseg; // normalized direction vector
2652 double ygrd = (yend - ylbend) / lseg;
2653 double zgrd = (zend - zlbend) / lseg;
2654 if (debugFn) {
2655 printf("\nelectron: %d\n", electron);
2656 printf("xlbend: %lg, ylbend: %lg, zlbend: %lg\n", xlbend, ylbend,
2657 zlbend);
2658 printf("xend: %lg, yend: %lg, zend: %lg\n", xend, yend, zend);
2659 printf("xgrd: %lg, ygrd: %lg, zgrd: %lg\n", xgrd, ygrd, zgrd);
2660 fprintf(ftmpEF, "#e: %d, label: %c, vol: %d\n", electron, label,
2661 vol);
2662 fprintf(ftmpEF, "%lg %lg %lg\n", xlbend, ylbend, zlbend);
2663 fprintf(ftmpEF, "%lg %lg %lg\n", xend, yend, zend);
2664 fprintf(ftmpEF, "#xgrd: %lg, ygrd: %lg, zgrd: %lg\n", xgrd, ygrd,
2665 zgrd);
2666 fprintf(ftmpEF, "\n");
2667 }
2668
2669 // determine which element gets this electron
2670 // At present, the logic is as follow:
2671 // Using the information on the last segment, find out which
2672 // primitive is pierced by it and at which point From the
2673 // intersection point, find out the element number on the primitive
2674 // in a local sense Using the information (start and end elements of
2675 // a given primitive) identify the element in a global sense. This
2676 // approach should be lot more efficient than checking intersection
2677 // element by element.
2678 // The intersection point is computed following algorithm
2679 // implemented in a matlab code (plane_imp_line_par_int_3d.m) Also
2680 // check which primitive in the list is the closet to the end point
2681
2682 double SumOfAngles;
2683 int PrimIntsctd = -1,
2684 EleIntsctd = -1; // intersected primitive & element
2685 int nearestprim = -1; // absurd value
2686 double dist = 1.0e6, mindist = 1.0e6; // absurdly high numbers
2687 // check all primitives
2688 for (int prim = 1; prim <= NbPrimitives; ++prim) {
2689 if (InterfaceType[prim] != 4)
2690 continue; // primitive not a dielectric
2691
2692 int intersect = 0, extrasect = 1; // worst of conditions
2693 int InPrim = 0, InEle = 0;
2694 if (debugFn)
2695 printf("prim: %d, mindist: %lg, nearestprim: %d\n", prim,
2696 mindist, nearestprim);
2697
2698 // Use two nodes at a time to get two independent vectors on
2699 // primitive Get cross-product of these two vector - normal to the
2700 // plane Note that the normal is already associated with the
2701 // primitive of type 3 and 4; 2 is wire and does not have any
2702 // associated normal
2703 if ((PrimType[prim] == 3) || (PrimType[prim] == 4)) {
2704 if (debugFn) {
2705 printf("prim: %d\n", prim);
2706 printf("vertex0: %lg, %lg, %lg\n", XVertex[prim][0],
2707 YVertex[prim][0], ZVertex[prim][0]);
2708 printf("vertex1: %lg, %lg, %lg\n", XVertex[prim][1],
2709 YVertex[prim][1], ZVertex[prim][1]);
2710 printf("vertex2: %lg, %lg, %lg\n", XVertex[prim][2],
2711 YVertex[prim][2], ZVertex[prim][2]);
2712 if (PrimType[prim] == 4) {
2713 printf("vertex3: %lg, %lg, %lg\n", XVertex[prim][3],
2714 YVertex[prim][3], ZVertex[prim][3]);
2715 }
2716 fprintf(ftmpEF, "#prim: %d\n", prim);
2717 fprintf(ftmpEF, "%lg %lg %lg\n", XVertex[prim][0],
2718 YVertex[prim][0], ZVertex[prim][0]);
2719 fprintf(ftmpEF, "%lg %lg %lg\n", XVertex[prim][1],
2720 YVertex[prim][1], ZVertex[prim][1]);
2721 fprintf(ftmpEF, "%lg %lg %lg\n", XVertex[prim][2],
2722 YVertex[prim][2], ZVertex[prim][2]);
2723 if (PrimType[prim] == 4) {
2724 fprintf(ftmpEF, "%lg %lg %lg\n", XVertex[prim][3],
2725 YVertex[prim][3], ZVertex[prim][3]);
2726 }
2727 fprintf(ftmpEF, "%lg %lg %lg\n", XVertex[prim][0],
2728 YVertex[prim][0], ZVertex[prim][0]);
2729 fprintf(ftmpEF, "\n");
2730 fflush(stdout);
2731 } // debugFn
2732
2733 // use a, b, c (normal is ai + bj + ck) at one of the nodes to
2734 // get d, ax + by + cz + d = 0 is the equation of the plane
2735 double a = XNorm[prim];
2736 double b = YNorm[prim];
2737 double c = ZNorm[prim];
2738 double d = -a * XVertex[prim][0] - b * YVertex[prim][0] -
2739 c * ZVertex[prim][0];
2740
2741 // distance of the end point to this primitve is
2742 dist = (xend * a + yend * b + zend * c + d) /
2743 sqrt(a * a + b * b + c * c);
2744 dist = fabs(dist); // if only magnitude is required
2745 if (prim == 1) {
2746 mindist = dist;
2747 nearestprim = prim;
2748 }
2749 if ((prim == 1) && debugFn)
2750 printf(
2751 "after prim == 1 check mindist: %lg, nearestprim: %d\n",
2752 mindist, nearestprim);
2753
2754 // Point of intersection
2755 // Algo as on p62 (pdf 81) of Vince - Geometry for Computer
2756 // Graphics 1.5.13 Intersection of a line and a plane Algorithm
2757 // as implemented in plne_imp_line_par_int_3d.m a (nx), b (ny),
2758 // c (nz), d are a, b, c, d vx, vy, vz are xgrd, ygrd and zgrd
2759 // tx, ty, tz are xlbend, ylbend, zlbend
2760 // In the present case, n and v are unit vectors
2761 double norm1 = sqrt(a * a + b * b + c * c);
2762 double norm2 = sqrt(xgrd * xgrd + ygrd * ygrd + zgrd * zgrd);
2763 double denom =
2764 a * xgrd + b * ygrd + c * zgrd; // (vec)n.(vec)v; if 0, ||
2765 double tol =
2766 1.0e-16; // CHECK: -8 in original code; sizes small here
2767 intersect = extrasect = 0;
2768
2769 if (debugFn) {
2770 printf("a, b, c, d, dist: %lg, %lg, %lg, %lg, %lg\n", a, b, c,
2771 d, dist);
2772 printf("vector n: ai + bj + ck\n");
2773 printf("vector v: xgrd, ygrd, zgrd: %lg, %lg, %lg\n", xgrd,
2774 ygrd, zgrd);
2775 printf("norm1, norm2, (vec n . vec v) denom: %lg, %lg, %lg\n",
2776 norm1, norm2, denom);
2777 printf("if vec n . vec v == 0, line and plane parallel\n");
2778 fflush(stdout);
2779 }
2780
2781 if (fabs(denom) < tol * norm1 * norm2) {
2782 // line parallel to the plane
2783 if (fabs(a * xlbend + b * ylbend + c * zlbend + d) <=
2784 1.0e-16) { // CHECK: was == 0.0 in original code
2785 intersect = 1;
2786 extrasect = 0; // line ends on the plane
2787 ptintsct.X = xlbend;
2788 ptintsct.Y = ylbend;
2789 ptintsct.Z = zlbend;
2790 } else {
2791 intersect = 0;
2792 extrasect = 1; // both wrong
2793 ptintsct.X = 0.0; // Wrong to assign 0 values
2794 ptintsct.Y =
2795 0.0; // However, they are never going to be used
2796 ptintsct.Z = 0.0; // since intersect is 0
2797 }
2798 if (debugFn) {
2799 printf("line and plane parallel ...\n");
2800 printf("intersect: %d, extrasect: %d\n", intersect,
2801 extrasect);
2802 printf("intersection point: %lg, %lg, %lg\n", ptintsct.X,
2803 ptintsct.Y, ptintsct.Z);
2804 } // if line and plane are parallel
2805 } else { // if they are not parallel, they must intersect
2806 intersect = 1;
2807 double t =
2808 -(a * xlbend + b * ylbend + c * zlbend + d) / denom;
2809
2810 // check whether t is less than the length of the segment
2811 // and in the correct direction
2812 // If not, then an extrapolated intersection is not of
2813 // interest
2814 if ((t < 0.0) ||
2815 (fabs(t) > fabs(lseg))) // wrong dirn or beyond end
2816 {
2817 extrasect = 1;
2818 ptintsct.X = xlbend + t * xgrd;
2819 ptintsct.Y = ylbend + t * ygrd;
2820 ptintsct.Z = zlbend + t * zgrd;
2821 } else {
2822 extrasect = 0;
2823 ptintsct.X = xlbend + t * xgrd;
2824 ptintsct.Y = ylbend + t * ygrd;
2825 ptintsct.Z = zlbend + t * zgrd;
2826 }
2827 if (debugFn) {
2828 printf("line and plane NOT parallel ...\n");
2829 printf("intersect: %d, extrasect: %d\n", intersect,
2830 extrasect);
2831 printf("intersection point: %lg, %lg, %lg\n", ptintsct.X,
2832 ptintsct.Y, ptintsct.Z);
2833 printf("t, lseg: %lg, %lg\n", t, lseg);
2834 printf(
2835 "for an interesting intersection, lseg > t > 0.0 "
2836 "...\n\n");
2837 fflush(stdout);
2838 } // must intersect
2839 } // if not parallel
2840 } // if PrimType is 3 or 4
2841 else { // this is a wire primitive - assume no charging up issues
2842 dist = -1.0; // an absurd negative distance
2843 intersect = 0;
2844 extrasect = 0;
2845 continue;
2846 } // else PrimType 3 or 4
2847
2848 if (dist < mindist) {
2849 mindist = dist;
2850 nearestprim = prim;
2851 }
2852 if (debugFn)
2853 printf("nearestprim: %d, mindist: %lg\n\n", nearestprim,
2854 mindist);
2855
2856 // implicit assumption: the first primitive that gets pierced by
2857 // the ray is the one that we want. There can be other primitives
2858 // that are pierced by the same ray. So, this logic should be
2859 // refined further
2860 if ((intersect == 1) && (extrasect == 0)) {
2861 // check whether the intersection point is within primitive
2862 // polygon
2863 int nvert = PrimType[prim];
2864 Point3D polynode[4];
2865 polynode[0].X = XVertex[prim][0];
2866 polynode[0].Y = YVertex[prim][0];
2867 polynode[0].Z = ZVertex[prim][0];
2868 polynode[1].X = XVertex[prim][1];
2869 polynode[1].Y = YVertex[prim][1];
2870 polynode[1].Z = ZVertex[prim][1];
2871 polynode[2].X = XVertex[prim][2];
2872 polynode[2].Y = YVertex[prim][2];
2873 polynode[2].Z = ZVertex[prim][2];
2874 if (PrimType[prim] == 4) {
2875 polynode[3].X = XVertex[prim][3];
2876 polynode[3].Y = YVertex[prim][3];
2877 polynode[3].Z = ZVertex[prim][3];
2878 }
2879 // printf("neBEMChkInPoly for primitive %d\n", prim);
2880 SumOfAngles = neBEMChkInPoly(nvert, polynode, ptintsct);
2881 if (fabs(fabs(SumOfAngles) - neBEMtwopi) <= 1.0e-8) {
2882 InPrim = 1;
2883 PrimIntsctd = prim;
2884 }
2885 if (debugFn) {
2886 // print polynode and InPrim
2887 printf("Prim: %d\n", prim);
2888 printf("ptintsct: %lg, %lg, %lg\n", ptintsct.X, ptintsct.Y,
2889 ptintsct.Z);
2890 printf("nvert: %d\n", nvert);
2891 printf("polynode0: %lg, %lg, %lg\n", polynode[0].X,
2892 polynode[0].Y, polynode[0].Z);
2893 printf("polynode1: %lg, %lg, %lg\n", polynode[1].X,
2894 polynode[1].Y, polynode[1].Z);
2895 printf("polynode2: %lg, %lg, %lg\n", polynode[2].X,
2896 polynode[2].Y, polynode[2].Z);
2897 if (nvert == 4) {
2898 printf("polynode3: %lg, %lg, %lg\n", polynode[3].X,
2899 polynode[3].Y, polynode[3].Z);
2900 }
2901 printf("SumOfAngles: %lg, InPrim: %d\n", SumOfAngles, InPrim);
2902 fflush(stdout);
2903 }
2904
2905 if (!InPrim && (prim != NbPrimitives)) {
2906 continue; // check next primitive
2907 }
2908
2909 // Once identified, check in which element belonging to this
2910 // primitive contains the point of intersection
2911 if (InPrim) {
2912 InEle = 0;
2913 for (int ele = ElementBgn[prim]; ele <= ElementEnd[prim];
2914 ++ele) {
2915 nvert = 0;
2916 if ((EleArr + ele - 1)->G.Type == 3) nvert = 3;
2917 if ((EleArr + ele - 1)->G.Type == 4) nvert = 4;
2918 if (!nvert) {
2920 "no vertex in element! ... neBEMKnownCharges ...\n");
2921 fclose(fPtEChUpMap);
2922 return -20;
2923 }
2924
2925 polynode[0].X = (EleArr + ele - 1)->G.Vertex[0].X;
2926 polynode[0].Y = (EleArr + ele - 1)->G.Vertex[0].Y;
2927 polynode[0].Z = (EleArr + ele - 1)->G.Vertex[0].Z;
2928 polynode[1].X = (EleArr + ele - 1)->G.Vertex[1].X;
2929 polynode[1].Y = (EleArr + ele - 1)->G.Vertex[1].Y;
2930 polynode[1].Z = (EleArr + ele - 1)->G.Vertex[1].Z;
2931 polynode[2].X = (EleArr + ele - 1)->G.Vertex[2].X;
2932 polynode[2].Y = (EleArr + ele - 1)->G.Vertex[2].Y;
2933 polynode[2].Z = (EleArr + ele - 1)->G.Vertex[2].Z;
2934 if (nvert == 4) {
2935 polynode[3].X = (EleArr + ele - 1)->G.Vertex[3].X;
2936 polynode[3].Y = (EleArr + ele - 1)->G.Vertex[3].Y;
2937 polynode[3].Z = (EleArr + ele - 1)->G.Vertex[3].Z;
2938 }
2939
2940 // printf("neBEMChkInPoly for element %d\n", ele);
2941 SumOfAngles = neBEMChkInPoly(nvert, polynode, ptintsct);
2942 if (fabs(fabs(SumOfAngles) - neBEMtwopi) <= 1.0e-8)
2943 InEle = 1;
2944 if (debugFn) {
2945 // print polynode and InEle
2946 printf("Ele: %d\n", ele);
2947 printf("ptintsct: %lg, %lg, %lg\n", ptintsct.X,
2948 ptintsct.Y, ptintsct.Z);
2949 printf("nvert: %d\n", nvert);
2950 printf("polynode0: %lg, %lg, %lg\n", polynode[0].X,
2951 polynode[0].Y, polynode[0].Z);
2952 printf("polynode1: %lg, %lg, %lg\n", polynode[1].X,
2953 polynode[1].Y, polynode[1].Z);
2954 printf("polynode2: %lg, %lg, %lg\n", polynode[2].X,
2955 polynode[2].Y, polynode[2].Z);
2956 if (nvert == 4) {
2957 printf("polynode3: %lg, %lg, %lg\n", polynode[3].X,
2958 polynode[3].Y, polynode[3].Z);
2959 }
2960 printf("SumOfAngles: %lg, InEle: %d\n", SumOfAngles,
2961 InEle);
2962 fflush(stdout);
2963 }
2964 if (InEle) {
2965 ptintsct.X = (EleArr + ele - 1)->G.Origin.X;
2966 ptintsct.Y = (EleArr + ele - 1)->G.Origin.Y;
2967 ptintsct.Z = (EleArr + ele - 1)->G.Origin.Z;
2968 // Associate this electron to the identified element
2969 EleIntsctd = ele;
2970 NbChUpEonEle[ele]++;
2971 fprintf(fPtEChUpMap, "%d %lg %lg %lg %d %d %d %d\n",
2972 electron, ptintsct.X, ptintsct.Y, ptintsct.Z,
2973 prim, InPrim, ele, InEle);
2974
2975 if (debugFn) {
2976 printf("# electron: %d\n", electron);
2977 printf("%lg %lg %lg\n", xlbend, ylbend, zlbend);
2978 printf("%lg %lg %lg\n", xend, yend, zend);
2979 printf("%lg, %lg, %lg\n", ptintsct.X, ptintsct.Y,
2980 ptintsct.Z);
2981 printf("# Associated primitive: %d\n", prim);
2982 printf(
2983 "# Associated element and origin: %d, %lg, %lg, "
2984 "%lg\n",
2985 ele, (EleArr + ele - 1)->G.Origin.X,
2986 (EleArr + ele - 1)->G.Origin.Y,
2987 (EleArr + ele - 1)->G.Origin.Z);
2988 printf("#NbChUpEonEle on element: %d\n",
2989 NbChUpEonEle[ele]);
2990 fprintf(ftmpEF, "#Element: %d\n", ele);
2991 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[0].X,
2992 polynode[0].Y, polynode[0].Z);
2993 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[1].X,
2994 polynode[1].Y, polynode[1].Z);
2995 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[2].X,
2996 polynode[2].Y, polynode[2].Z);
2997 if (nvert == 4) {
2998 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[3].X,
2999 polynode[3].Y, polynode[3].Z);
3000 }
3001 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[0].X,
3002 polynode[0].Y, polynode[0].Z);
3003 fprintf(ftmpEF, "\n");
3004 fflush(stdout);
3005 }
3006 break; // desired element has been found!
3007 }
3008 } // for all elements on this primitive
3009
3010 if (InEle)
3011 break;
3012 else {
3014 "Element not identified ... neBEMKnownCharges\n");
3015 return -2;
3016 }
3017 } // if InPrim
3018 } // if intersection and no extrasection
3019
3020 if ((InPrim) && (intersect) && (!extrasect) && (InEle)) {
3021 // all satisfied
3022 // do not check any further primitive for this electron
3023 break;
3024 }
3025
3026 // If, after checking all the primitives, no interstion is found
3027 // valid
3028 if (prim ==
3029 (NbPrimitives)) // end of the list and no intersection
3030 {
3031 int nvert;
3032 Point3D polynode[4];
3033 int nearestele = ElementBgn[nearestprim];
3034 double distele = 1.0e6;
3035 double mindistele = 1.0e6; // absurdly high value
3036
3037 if (debugFn) {
3038 printf("prim == (NbPrimitives) ... checking nearest ...\n");
3039 printf("nearestprim: %d, mindist: %lg\n", nearestprim,
3040 mindist);
3041 }
3042
3043 if (mindist <= 10.0e-6) {
3044 PrimIntsctd = nearestprim;
3045 InPrim = 1;
3046 } else {
3047 InPrim = 0;
3048 InEle = 0;
3049 break;
3050 }
3051 // check all elements
3052 for (int ele = ElementBgn[nearestprim];
3053 ele <= ElementEnd[nearestprim]; ++ele) {
3054 nvert = 0;
3055 if ((EleArr + ele - 1)->G.Type == 3) nvert = 3;
3056 if ((EleArr + ele - 1)->G.Type == 4) nvert = 4;
3057 if (!nvert) {
3059 "no vertex element! ... neBEMKnownCharges ...\n");
3060 return -20;
3061 }
3062
3063 /*
3064 polynode[0].X = (EleArr+ele-1)->G.Vertex[0].X;
3065 polynode[0].Y = (EleArr+ele-1)->G.Vertex[0].Y;
3066 polynode[0].Z = (EleArr+ele-1)->G.Vertex[0].Z;
3067 polynode[1].X = (EleArr+ele-1)->G.Vertex[1].X;
3068 polynode[1].Y = (EleArr+ele-1)->G.Vertex[1].Y;
3069 polynode[1].Z = (EleArr+ele-1)->G.Vertex[1].Z;
3070 polynode[2].X = (EleArr+ele-1)->G.Vertex[2].X;
3071 polynode[2].Y = (EleArr+ele-1)->G.Vertex[2].Y;
3072 polynode[2].Z = (EleArr+ele-1)->G.Vertex[2].Z;
3073 if (nvert == 4) {
3074 polynode[3].X = (EleArr+ele-1)->G.Vertex[3].X;
3075 polynode[3].Y = (EleArr+ele-1)->G.Vertex[3].Y;
3076 polynode[3].Z = (EleArr+ele-1)->G.Vertex[3].Z;
3077 }
3078 Vector3D v01, v12, elenorm, unitelenorm;
3079 v01.X = polynode[1].X - polynode[0].X;
3080 v01.Y = polynode[1].Y - polynode[0].Y;
3081 v01.Z = polynode[1].Z - polynode[0].Z;
3082 v12.X = polynode[2].X - polynode[1].X;
3083 v12.Y = polynode[2].Y - polynode[1].Y;
3084 v12.Z = polynode[2].Z - polynode[1].Z;
3085 elenorm = Vector3DCrossProduct(&v01, &v12);
3086 unitelenorm = UnitVector3D(&elenorm);
3087
3088 if ((nvert == 3) || (nvert == 4)) {
3089 if (debugFn) {
3090 printf("nearestprim: %d, element: %d\n",
3091 nearestprim, ele);
3092 printf("vertex0: %lg, %lg, %lg\n", polynode[0].X,
3093 polynode[0].Y, polynode[0].Z);
3094 printf("vertex1: %lg, %lg, %lg\n", polynode[1].X,
3095 polynode[1].Y, polynode[1].Z);
3096 printf("vertex2: %lg, %lg, %lg\n", polynode[2].X,
3097 polynode[2].Y, polynode[2].Z);
3098 if (PrimType[prim] == 4) {
3099 printf("vertex3: %lg, %lg, %lg\n", polynode[3].X,
3100 polynode[3].Y, polynode[3].Z);
3101 }
3102 fprintf(ftmpEF, "#nearestprim: %d, element: %d\n",
3103 nearestprim, ele);
3104 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[0].X,
3105 polynode[0].Y, polynode[0].Z);
3106 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[1].X,
3107 polynode[1].Y, polynode[1].Z);
3108 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[2].X,
3109 polynode[2].Y, polynode[2].Z);
3110 if (PrimType[prim] == 4) {
3111 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[3].X,
3112 polynode[3].Y, polynode[3].Z);
3113 }
3114 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[0].X,
3115 polynode[0].Y, polynode[0].Z);
3116 fprintf(ftmpEF, "\n");
3117 fflush(stdout);
3118 } // debugFn
3119 // use a, b, c (normal is ai + bj + ck)
3120 // at one of the nodes to get d
3121 // ax + by + cz + d = 0 is the equation of the plane
3122 double a = unitelenorm.X;
3123 double b = unitelenorm.Y;
3124 double c = unitelenorm.Z;
3125 double d = - a * polynode[0].X - b * polynode[0].Y - c * polynode[0].Z;
3126 // distance of the end point to this primitve is
3127 distele = (xend * a + yend * b + zend * c + d) /
3128 sqrt(a * a + b * b + c * c);
3129 distele = fabs(distele); // if only magnitude is required
3130 */
3131
3132 Vector3D eleOrigin;
3133 eleOrigin.X = (EleArr + ele - 1)->G.Origin.X;
3134 eleOrigin.Y = (EleArr + ele - 1)->G.Origin.Y;
3135 eleOrigin.Z = (EleArr + ele - 1)->G.Origin.Z;
3136 distele = (eleOrigin.X - xend) * (eleOrigin.X - xend) +
3137 (eleOrigin.Y - yend) * (eleOrigin.Y - yend) +
3138 (eleOrigin.Z - zend) * (eleOrigin.Z - zend);
3139 distele = pow(distele, 0.5);
3140
3141 if (ele == ElementBgn[nearestprim]) {
3142 mindistele = distele;
3143 nearestele = ele;
3144 }
3145 if (distele < mindistele) {
3146 mindistele = distele;
3147 nearestele = ele;
3148 }
3149
3150 if (debugFn) {
3151 // printf("a, b, c, d, dist: %lg, %lg, %lg, %lg, %lg\n",
3152 // a, b, c, d, dist);
3153 // printf("vector n: ai + bj + ck\n");
3154 // printf("vector v: xgrd, ygrd, zgrd: %lg, %lg, %lg\n",
3155 // xgrd, ygrd, zgrd);
3156 printf(
3157 "distele: %lg, mindistele: %lg,from nearest ele "
3158 "origin: %d\n",
3159 distele, mindistele, nearestele);
3160 fflush(stdout);
3161 }
3162
3163 // } // if PrimType is 3 or 4
3164 } // for elements in nearestprim
3165
3166 // if(mindistele <= 10.0e-6)
3167 // {
3168 EleIntsctd = nearestele;
3169 InEle = 1;
3170 ptintsct.X = (EleArr + EleIntsctd - 1)->G.Origin.X;
3171 ptintsct.Y = (EleArr + EleIntsctd - 1)->G.Origin.Y;
3172 ptintsct.Z = (EleArr + EleIntsctd - 1)->G.Origin.Z;
3173 NbChUpEonEle[EleIntsctd]++;
3174
3175 fprintf(fPtEChUpMap, "%d %lg %lg %lg %d %d %d %d\n", electron,
3176 ptintsct.X, ptintsct.Y, ptintsct.Z, PrimIntsctd, InPrim,
3177 EleIntsctd, InEle);
3178 // } // if mindistele
3179
3180 if (debugFn) {
3181 printf("# electron: %d\n", electron);
3182 printf("%lg %lg %lg\n", xlbend, ylbend, zlbend);
3183 printf("%lg %lg %lg\n", xend, yend, zend);
3184 printf("%lg, %lg, %lg\n", ptintsct.X, ptintsct.Y, ptintsct.Z);
3185 printf("# Associated primitive: %d\n", PrimIntsctd);
3186 printf("# Associated element and origin: %d, %lg, %lg, %lg\n",
3187 EleIntsctd, (EleArr + EleIntsctd - 1)->G.Origin.X,
3188 (EleArr + EleIntsctd - 1)->G.Origin.Y,
3189 (EleArr + EleIntsctd - 1)->G.Origin.Z);
3190 printf("#NbChUpEonEle on element: %d\n",
3191 NbChUpEonEle[EleIntsctd]);
3192 fflush(stdout);
3193
3194 fprintf(ftmpEF, "#Element: %d\n", EleIntsctd);
3195 polynode[0].X = (EleArr + EleIntsctd - 1)->G.Vertex[0].X;
3196 polynode[0].Y = (EleArr + EleIntsctd - 1)->G.Vertex[0].Y;
3197 polynode[0].Z = (EleArr + EleIntsctd - 1)->G.Vertex[0].Z;
3198 polynode[1].X = (EleArr + EleIntsctd - 1)->G.Vertex[1].X;
3199 polynode[1].Y = (EleArr + EleIntsctd - 1)->G.Vertex[1].Y;
3200 polynode[1].Z = (EleArr + EleIntsctd - 1)->G.Vertex[1].Z;
3201 polynode[2].X = (EleArr + EleIntsctd - 1)->G.Vertex[2].X;
3202 polynode[2].Y = (EleArr + EleIntsctd - 1)->G.Vertex[2].Y;
3203 polynode[2].Z = (EleArr + EleIntsctd - 1)->G.Vertex[2].Z;
3204 if (nvert == 4) {
3205 polynode[3].X = (EleArr + EleIntsctd - 1)->G.Vertex[3].X;
3206 polynode[3].Y = (EleArr + EleIntsctd - 1)->G.Vertex[3].Y;
3207 polynode[3].Z = (EleArr + EleIntsctd - 1)->G.Vertex[3].Z;
3208 }
3209 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[0].X, polynode[0].Y,
3210 polynode[0].Z);
3211 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[1].X, polynode[1].Y,
3212 polynode[1].Z);
3213 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[2].X, polynode[2].Y,
3214 polynode[2].Z);
3215 if (nvert == 4) {
3216 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[3].X,
3217 polynode[3].Y, polynode[3].Z);
3218 }
3219 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[0].X, polynode[0].Y,
3220 polynode[0].Z);
3221 fprintf(ftmpEF, "\n");
3222 } // debug
3223 } // if prim == NbPrimitives
3224
3225 } // for all primitives // just not those on the volume
3226
3227 if (debugFn)
3228 printf("writing file for checking electron positions\n");
3229
3230 if (debugFn) // check electron positions, volume primitives and
3231 // elements
3232 {
3233 char elecposdbg[256], enbstr[10];
3234 sprintf(enbstr, "%d", electron);
3235 strcpy(elecposdbg, "/tmp/Electron");
3236 strcat(elecposdbg, enbstr);
3237 strcat(elecposdbg, ".out");
3238 FILE *fepd = fopen(elecposdbg, "w");
3239 if (fepd == NULL) {
3240 printf(
3241 "cannot open writable file to debug electron positions "
3242 "...\n");
3243 printf("returning ...\n");
3244 return -111;
3245 }
3246 // write electron number, end, lbend, volume, primitive, elements,
3247 // intxn
3248 fprintf(fepd, "#electron: %d %d\n", enb,
3249 electron); // should print same
3250 fprintf(fepd, "#last but end position:\n");
3251 fprintf(fepd, "%lg %lg %lg\n", xlbend, ylbend, zlbend);
3252 fprintf(fepd, "#end position:\n");
3253 fprintf(fepd, "%lg %lg %lg\n\n", xend, yend, zend);
3254 fprintf(fepd, "#intersected primitive number: %d\n", PrimIntsctd);
3255 if (PrimIntsctd >= 1) {
3256 fprintf(fepd, "#PrimType: %d\n", PrimType[PrimIntsctd]);
3257 fprintf(fepd, "#prim vertices:\n");
3258 fprintf(fepd, "%lg %lg %lg\n", XVertex[PrimIntsctd][0],
3259 YVertex[PrimIntsctd][0], ZVertex[PrimIntsctd][0]);
3260 fprintf(fepd, "%lg %lg %lg\n", XVertex[PrimIntsctd][1],
3261 YVertex[PrimIntsctd][1], ZVertex[PrimIntsctd][1]);
3262 fprintf(fepd, "%lg %lg %lg\n", XVertex[PrimIntsctd][2],
3263 YVertex[PrimIntsctd][2], ZVertex[PrimIntsctd][2]);
3264 if (PrimType[PrimIntsctd] == 4) {
3265 fprintf(fepd, "%lg %lg %lg\n", XVertex[PrimIntsctd][3],
3266 YVertex[PrimIntsctd][3], ZVertex[PrimIntsctd][3]);
3267 }
3268 fprintf(fepd, "%lg %lg %lg\n", XVertex[PrimIntsctd][0],
3269 YVertex[PrimIntsctd][0], ZVertex[PrimIntsctd][0]);
3270 fprintf(fepd, "\n");
3271
3272 fprintf(fepd, "#ptintsct:\n");
3273 fprintf(fepd, "%lg %lg %lg\n", ptintsct.X, ptintsct.Y,
3274 ptintsct.Z);
3275 fprintf(fepd, "\n");
3276 }
3277 fprintf(fepd, "#intersected element number: %d\n", EleIntsctd);
3278 if (EleIntsctd >= 1) {
3279 int gtype = (EleArr + EleIntsctd - 1)->G.Type;
3280 fprintf(fepd, "#EleType: %d\n", gtype);
3281 fprintf(fepd, "#element vertices:\n");
3282 double x0 = (EleArr + EleIntsctd - 1)->G.Vertex[0].X;
3283 double y0 = (EleArr + EleIntsctd - 1)->G.Vertex[0].Y;
3284 double z0 = (EleArr + EleIntsctd - 1)->G.Vertex[0].Z;
3285 double x1 = (EleArr + EleIntsctd - 1)->G.Vertex[1].X;
3286 double y1 = (EleArr + EleIntsctd - 1)->G.Vertex[1].Y;
3287 double z1 = (EleArr + EleIntsctd - 1)->G.Vertex[1].Z;
3288 double x2 = (EleArr + EleIntsctd - 1)->G.Vertex[2].X;
3289 double y2 = (EleArr + EleIntsctd - 1)->G.Vertex[2].Y;
3290 double z2 = (EleArr + EleIntsctd - 1)->G.Vertex[2].Z;
3291 fprintf(fepd, "%lg %lg %lg\n", x0, y0, z0);
3292 fprintf(fepd, "%lg %lg %lg\n", x1, y1, z1);
3293 fprintf(fepd, "%lg %lg %lg\n", x2, y2, z2);
3294 if (gtype == 4) {
3295 double x3 = (EleArr + EleIntsctd - 1)->G.Vertex[3].X;
3296 double y3 = (EleArr + EleIntsctd - 1)->G.Vertex[3].Y;
3297 double z3 = (EleArr + EleIntsctd - 1)->G.Vertex[3].Z;
3298 fprintf(fepd, "%lg %lg %lg\n", x3, y3, z3);
3299 }
3300 fprintf(fepd, "%lg %lg %lg\n", x0, y0, z0);
3301 fprintf(fepd, "\n");
3302
3303 fprintf(fepd, "#ptintsct:\n");
3304 fprintf(fepd, "%lg %lg %lg\n", ptintsct.X, ptintsct.Y,
3305 ptintsct.Z);
3306 fprintf(fepd, "\n");
3307 }
3308
3309 fclose(fepd);
3310 } // if 1
3311 if (debugFn)
3312 printf("done writing file for checking electron positions\n");
3313 } // for all the electrons
3314 fclose(fPtEChUpMap);
3315 if (debugFn) printf("done for all electrons\n");
3316
3317 FILE *fEleEChUpMap = fopen("EleEChUpMap.out", "w");
3318 if (fEleEChUpMap == NULL) {
3319 printf("cannot open EleEChUpMap.out file for writing ...\n");
3320 return 111;
3321 }
3322 for (int ele = 1; ele <= NbElements; ++ele) {
3323 (EleArr + ele - 1)->Assigned +=
3324 ChUpFactor * Q_E * NbChUpEonEle[ele] / (EleArr + ele - 1)->G.dA;
3325 fprintf(fEleEChUpMap, "%d %lg %lg %lg %d %lg\n", ele,
3326 (EleArr + ele - 1)->G.Origin.X,
3327 (EleArr + ele - 1)->G.Origin.Y,
3328 (EleArr + ele - 1)->G.Origin.Z, NbChUpEonEle[ele],
3329 (EleArr + ele - 1)->Assigned);
3330 }
3331 fclose(fEleEChUpMap);
3332
3333 fclose(ftmpEF);
3334 free(NbChUpEonEle);
3335 } // Calculation for electrons ends
3336
3337 { // Calculation for ions
3338 FILE *fptrChargingUpIFile = fopen(ChargingUpIFile, "r");
3339 if (fptrChargingUpIFile == NULL) {
3340 neBEMMessage("ChargingUpI file absent ... returning\n");
3341 return -10;
3342 }
3343 int NbOfI = neBEMGetNbOfLines(ChargingUpIFile);
3344 if (NbOfI <= 1) {
3345 neBEMMessage("Too few lines in ChargingUpI ... returning\n");
3346 fclose(fptrChargingUpIFile);
3347 return -11;
3348 }
3349 // initialize
3350 NbChUpIonEle = (int *)malloc((NbElements + 1) * sizeof(int));
3351 if (NbChUpIonEle == NULL) {
3352 neBEMMessage("Memory allocation failed ... returning\n");
3353 fclose(fptrChargingUpIFile);
3354 return -11;
3355 }
3356 for (int ele = 0; ele <= NbElements; ++ele) {
3357 // CHECK!!! ele limit starts from 0 but all other from 1 to ...
3358 NbChUpIonEle[ele] = 0;
3359 }
3360
3361 // read the header line
3362 char header[256];
3363 fgets(header, 256, fptrChargingUpIFile);
3364
3365 --NbOfI; // one line was for the header
3366 if (debugFn) printf("No. of ions: %d\n", NbOfI);
3367 const char *tmpIFile = "/tmp/IonTempFile.out";
3368 FILE *ftmpIF = fopen(tmpIFile, "w");
3369 if (ftmpIF == NULL) {
3370 printf("cannot open temporary ion output file ... returning ...\n");
3371 free(NbChUpEonEle);
3372 return -100;
3373 }
3374 FILE *fPtIChUpMap = fopen("PtIChUpMap.out", "w");
3375 if (fPtIChUpMap == NULL) {
3376 printf("cannot open PtIChUpMap.out file for writing ...\n");
3377 fclose(ftmpIF);
3378 free(NbChUpEonEle);
3379 return 110;
3380 }
3381
3382 char label;
3383 int inb, vol; // label, volume and ion number
3384 double xlbend, ylbend, zlbend; // lbend == Last But END
3385 double xend, yend, zend;
3386 Point3D ptintsct; // each ion likely to have an intersection point
3387 for (int ion = 1; ion <= NbOfI; ++ion) {
3388 fscanf(fptrChargingUpIFile, "%c %d %d %lg %lg %lg %lg %lg %lg\n",
3389 &label, &vol, &inb, &xlbend, &xend, &ylbend, &yend, &zlbend,
3390 &zend);
3391 xlbend *= 0.01;
3392 xend *= 0.01;
3393 ylbend *= 0.01;
3394 yend *= 0.01;
3395 zlbend *= 0.01;
3396 zend *= 0.01;
3397 ptintsct.X = 0.0;
3398 ptintsct.Y = 0.0;
3399 ptintsct.Z = 0.0; // initialize
3400
3401 // find the parametric equation of this last segment
3402 // if xend < xlbend, swap the directions
3403 // This has not been mentioned as mandatory in Vince's book
3404 // "Geometry for Computer Graphics", but implied in the book "A
3405 // Programmer's Geometry"
3406 if (xend < xlbend) // not implemented now
3407 {
3408 }
3409 double lseg = (xend - xlbend) * (xend - xlbend) +
3410 (yend - ylbend) * (yend - ylbend) +
3411 (zend - zlbend) * (zend - zlbend);
3412 lseg = sqrt(lseg);
3413 double xgrd =
3414 (xend - xlbend) / lseg; // normalized direction vector
3415 double ygrd = (yend - ylbend) / lseg;
3416 double zgrd = (zend - zlbend) / lseg;
3417 if (debugFn) {
3418 printf("\nion: %d\n", ion);
3419 printf("xlbend: %lg, ylbend: %lg, zlbend: %lg\n", xlbend, ylbend,
3420 zlbend);
3421 printf("xend: %lg, yend: %lg, zend: %lg\n", xend, yend, zend);
3422 printf("xgrd: %lg, ygrd: %lg, zgrd: %lg\n", xgrd, ygrd, zgrd);
3423 fprintf(ftmpIF, "#e: %d, label: %c, vol: %d\n", ion, label, vol);
3424 fprintf(ftmpIF, "%lg %lg %lg\n", xlbend, ylbend, zlbend);
3425 fprintf(ftmpIF, "%lg %lg %lg\n", xend, yend, zend);
3426 fprintf(ftmpIF, "#xgrd: %lg, ygrd: %lg, zgrd: %lg\n", xgrd, ygrd,
3427 zgrd);
3428 fprintf(ftmpIF, "\n");
3429 }
3430
3431 // determine which element gets this electron
3432 // At present, the logic is as follow:
3433 // Using the information on the last segment, find out which
3434 // primitive is pierced by it and at which point From the
3435 // intersection point, find out the element number on the primitive
3436 // in a local sense Using the information (start and end elements of
3437 // a given primitive) identify the element in a global sense. This
3438 // approach should be lot more efficient than checking intersection
3439 // element by element.
3440 // The intersection point is computed following algorithm
3441 // implemented in a matlab code (plane_imp_line_par_int_3d.m) Also
3442 // check which primitive in the list is the closet to the end point
3443
3444 int PrimIntsctd = -1,
3445 EleIntsctd = -1; // intersected primitive & element
3446 int nearestprim = -1; // absurd value
3447 double dist = 1.0e6, mindist = 1.0e6; // absurdly high numbers
3448 double SumOfAngles;
3449 // check all primitives
3450 for (int prim = 1; prim <= NbPrimitives; ++prim) {
3451 if (InterfaceType[prim] != 4)
3452 continue; // primitive not a dielectric
3453
3454 int intersect = 0, extrasect = 1; // worst of conditions
3455 int InPrim = 0, InEle = 0;
3456 if (debugFn)
3457 printf("prim: %d, mindist: %lg, nearestprim: %d\n", prim,
3458 mindist, nearestprim);
3459
3460 // get the primitive nodes
3461
3462 // Use two nodes at a time to get two independent vectors on
3463 // primitive Get cross-product of these two vector - normal to the
3464 // plane Note that the normal is already associated with the
3465 // primitive of type 3 and 4; 2 is wire and does not have any
3466 // associated normal
3467 if ((PrimType[prim] == 3) || (PrimType[prim] == 4)) {
3468 if (debugFn) {
3469 printf("prim: %d\n", prim);
3470 printf("vertex0: %lg, %lg, %lg\n", XVertex[prim][0],
3471 YVertex[prim][0], ZVertex[prim][0]);
3472 printf("vertex1: %lg, %lg, %lg\n", XVertex[prim][1],
3473 YVertex[prim][1], ZVertex[prim][1]);
3474 printf("vertex2: %lg, %lg, %lg\n", XVertex[prim][2],
3475 YVertex[prim][2], ZVertex[prim][2]);
3476 if (PrimType[prim] == 4) {
3477 printf("vertex3: %lg, %lg, %lg\n", XVertex[prim][3],
3478 YVertex[prim][3], ZVertex[prim][3]);
3479 }
3480 fprintf(ftmpIF, "#prim: %d\n", prim);
3481 fprintf(ftmpIF, "%lg %lg %lg\n", XVertex[prim][0],
3482 YVertex[prim][0], ZVertex[prim][0]);
3483 fprintf(ftmpIF, "%lg %lg %lg\n", XVertex[prim][1],
3484 YVertex[prim][1], ZVertex[prim][1]);
3485 fprintf(ftmpIF, "%lg %lg %lg\n", XVertex[prim][2],
3486 YVertex[prim][2], ZVertex[prim][2]);
3487 if (PrimType[prim] == 4) {
3488 fprintf(ftmpIF, "%lg %lg %lg\n", XVertex[prim][3],
3489 YVertex[prim][3], ZVertex[prim][3]);
3490 }
3491 fprintf(ftmpIF, "%lg %lg %lg\n", XVertex[prim][0],
3492 YVertex[prim][0], ZVertex[prim][0]);
3493 fprintf(ftmpIF, "\n");
3494 fflush(stdout);
3495 } // debugFn
3496
3497 // use a, b, c (normal is ai + bj + ck) at one of the nodes to
3498 // get d ax + by + cz + d = 0 is the equation of the plane
3499 double d = -XNorm[prim] * XVertex[prim][0] -
3500 YNorm[prim] * YVertex[prim][0] -
3501 ZNorm[prim] * ZVertex[prim][0];
3502
3503 // distance of the end point to this primitve is
3504 dist =
3505 (xend * XNorm[prim] + yend * YNorm[prim] +
3506 zend * ZNorm[prim] + d) /
3507 sqrt(XNorm[prim] * XNorm[prim] + YNorm[prim] * YNorm[prim] +
3508 ZNorm[prim] * ZNorm[prim]);
3509 dist = fabs(dist); // if only magnitude is required
3510 if (prim == 1) {
3511 mindist = dist;
3512 nearestprim = prim;
3513 }
3514 if ((prim == 1) && debugFn)
3515 printf(
3516 "after prim == 1 check mindist: %lg, nearestprim: %d\n",
3517 mindist, nearestprim);
3518
3519 // Point of intersection
3520 // Algo as on p62 (pdf 81) of Vince - Geometry for Computer
3521 // Graphics 1.5.13 Intersection of a line and a plane Algorithm
3522 // as implemented in plne_imp_line_par_int_3d.m a (nx), b (ny),
3523 // c (nz), d are a, b, c, d vx, vy, vz are xgrd, ygrd and zgrd
3524 // tx, ty, tz are xlbend, ylbend, zlbend
3525 // In the present case, n and v are unit vectors
3526 double a = XNorm[prim];
3527 double b = YNorm[prim];
3528 double c = ZNorm[prim];
3529 double norm1 = sqrt(a * a + b * b + c * c);
3530 double norm2 = sqrt(xgrd * xgrd + ygrd * ygrd + zgrd * zgrd);
3531 double denom =
3532 a * xgrd + b * ygrd + c * zgrd; // (vec)n.(vec)v; if 0, ||
3533 double tol =
3534 1.0e-12; // CHECK: -8 in original code; sizes small here
3535 intersect = extrasect = 0;
3536
3537 if (debugFn) {
3538 printf("a, b, c, d, dist: %lg, %lg, %lg, %lg, %lg\n", a, b, c,
3539 d, dist);
3540 printf("vector n: ai + bj + ck\n");
3541 printf("vector v: xgrd, ygrd, zgrd: %lg, %lg, %lg\n", xgrd,
3542 ygrd, zgrd);
3543 printf("norm1, norm2, (vec n . vec v) denom: %lg, %lg, %lg\n",
3544 norm1, norm2, denom);
3545 printf("if vec n . vec v == 0, line and plane parallel\n");
3546 fflush(stdout);
3547 }
3548
3549 if (fabs(denom) < tol * norm1 * norm2) {
3550 // line parallel to the plane
3551 if (a * xlbend + b * ylbend + c * zlbend + d ==
3552 0.0) // CHECK == for float
3553 {
3554 intersect = 1;
3555 extrasect = 0;
3556 ptintsct.X = xlbend;
3557 ptintsct.Y = ylbend;
3558 ptintsct.Z = zlbend;
3559 } else {
3560 intersect = 0;
3561 extrasect = 0;
3562 ptintsct.X = 0.0; // Wrong to assign 0 values
3563 ptintsct.Y =
3564 0.0; // However, they are never going to be used
3565 ptintsct.Z = 0.0; // since intersect is 0
3566 }
3567 if (debugFn) {
3568 printf("line and plane parallel ...\n");
3569 printf("intersect: %d, extrasect: %d\n", intersect,
3570 extrasect);
3571 printf("intersection point: %lg, %lg, %lg\n", ptintsct.X,
3572 ptintsct.Y, ptintsct.Z);
3573 } // if line and plane are parallel
3574 } else { // if they are not parallel, they must intersect
3575 intersect = 1;
3576 double t =
3577 -(a * xlbend + b * ylbend + c * zlbend + d) / denom;
3578
3579 // check whether t is less than the length of the segment
3580 // and in the correct direction
3581 // If not, then an extrapolated intersection is not of
3582 // interest
3583 if ((t < 0.0) || (fabs(t) > fabs(lseg))) {
3584 extrasect = 1;
3585 ptintsct.X = xlbend + t * xgrd;
3586 ptintsct.Y = ylbend + t * ygrd;
3587 ptintsct.Z = zlbend + t * zgrd;
3588 } else {
3589 extrasect = 0;
3590 ptintsct.X = xlbend + t * xgrd;
3591 ptintsct.Y = ylbend + t * ygrd;
3592 ptintsct.Z = zlbend + t * zgrd;
3593 }
3594 if (debugFn) {
3595 printf("line and plane NOT parallel ...\n");
3596 printf("intersect: %d, extrasect: %d\n", intersect,
3597 extrasect);
3598 printf("intersection point: %lg, %lg, %lg\n", ptintsct.X,
3599 ptintsct.Y, ptintsct.Z);
3600 printf("t, lseg: %lg, %lg\n", t, lseg);
3601 printf(
3602 "for an interesting intersection, lseg > t > 0.0 "
3603 "...\n\n");
3604 fflush(stdout);
3605 } // must intersect
3606 } // if not parallel
3607 } // if PrimType is 3 or 4
3608 else { // this is a wire primitive - assume no charging up issues
3609 dist = -1.0; // an absurd negative distance
3610 intersect = 0;
3611 extrasect = 0;
3612 } // else PrimType 3 or 4
3613
3614 if (dist < mindist) {
3615 mindist = dist;
3616 nearestprim = prim;
3617 }
3618
3619 // implicit assumption: the first primitive that gets pierced by
3620 // the ray is the one that we want. There can be other primitives
3621 // that are pierced by the same ray. So, this logic should be
3622 // refined further
3623 if ((intersect == 1) && (extrasect == 0)) {
3624 // check whether the intersection point is within primitive
3625 // polygon
3626 int nvert = PrimType[prim];
3627 Point3D polynode[4];
3628 polynode[0].X = XVertex[prim][0];
3629 polynode[0].Y = YVertex[prim][0];
3630 polynode[0].Z = ZVertex[prim][0];
3631 polynode[1].X = XVertex[prim][1];
3632 polynode[1].Y = YVertex[prim][1];
3633 polynode[1].Z = ZVertex[prim][1];
3634 polynode[2].X = XVertex[prim][2];
3635 polynode[2].Y = YVertex[prim][2];
3636 polynode[2].Z = ZVertex[prim][2];
3637 if (PrimType[prim] == 4) {
3638 polynode[3].X = XVertex[prim][3];
3639 polynode[3].Y = YVertex[prim][3];
3640 polynode[3].Z = ZVertex[prim][3];
3641 }
3642 // printf("neBEMChkInPoly for primitive %d\n", prim);
3643 SumOfAngles = neBEMChkInPoly(nvert, polynode, ptintsct);
3644 if (fabs(fabs(SumOfAngles) - neBEMtwopi) <= 1.0e-8) {
3645 InPrim = 1;
3646 PrimIntsctd = prim;
3647 }
3648 if (debugFn) {
3649 // print polynode and InPrim
3650 printf("Prim: %d\n", prim);
3651 printf("ptintsct: %lg, %lg, %lg\n", ptintsct.X, ptintsct.Y,
3652 ptintsct.Z);
3653 printf("nvert: %d\n", nvert);
3654 printf("polynode0: %lg, %lg, %lg\n", polynode[0].X,
3655 polynode[0].Y, polynode[0].Z);
3656 printf("polynode1: %lg, %lg, %lg\n", polynode[1].X,
3657 polynode[1].Y, polynode[1].Z);
3658 printf("polynode2: %lg, %lg, %lg\n", polynode[2].X,
3659 polynode[2].Y, polynode[2].Z);
3660 if (nvert == 4) {
3661 printf("polynode3: %lg, %lg, %lg\n", polynode[3].X,
3662 polynode[3].Y, polynode[3].Z);
3663 }
3664 printf("SumOfAngles: %lg, InPrim: %d\n", SumOfAngles, InPrim);
3665 fflush(stdout);
3666 }
3667 if (!InPrim) continue; // check next primitive
3668
3669 // Once identified, check in which element belonging to this
3670 // primitive contains the point of intersection
3671 InEle = 0;
3672 for (int ele = ElementBgn[prim]; ele <= ElementEnd[prim];
3673 ++ele) {
3674 nvert = 0;
3675 if ((EleArr + ele - 1)->G.Type == 3) nvert = 3;
3676 if ((EleArr + ele - 1)->G.Type == 4) nvert = 4;
3677 if (!nvert) {
3679 "no vertex in element! ... neBEMKnownCharges ...\n");
3680 if (fPtIChUpMap) fclose(fPtIChUpMap);
3681 return -20;
3682 }
3683
3684 polynode[0].X = (EleArr + ele - 1)->G.Vertex[0].X;
3685 polynode[0].Y = (EleArr + ele - 1)->G.Vertex[0].Y;
3686 polynode[0].Z = (EleArr + ele - 1)->G.Vertex[0].Z;
3687 polynode[1].X = (EleArr + ele - 1)->G.Vertex[1].X;
3688 polynode[1].Y = (EleArr + ele - 1)->G.Vertex[1].Y;
3689 polynode[1].Z = (EleArr + ele - 1)->G.Vertex[1].Z;
3690 polynode[2].X = (EleArr + ele - 1)->G.Vertex[2].X;
3691 polynode[2].Y = (EleArr + ele - 1)->G.Vertex[2].Y;
3692 polynode[2].Z = (EleArr + ele - 1)->G.Vertex[2].Z;
3693 if (nvert == 4) {
3694 polynode[3].X = (EleArr + ele - 1)->G.Vertex[3].X;
3695 polynode[3].Y = (EleArr + ele - 1)->G.Vertex[3].Y;
3696 polynode[3].Z = (EleArr + ele - 1)->G.Vertex[3].Z;
3697 }
3698
3699 // printf("neBEMChkInPoly for element %d\n", ele);
3700 SumOfAngles = neBEMChkInPoly(nvert, polynode, ptintsct);
3701 if (fabs(fabs(SumOfAngles) - neBEMtwopi) <= 1.0e-8) InEle = 1;
3702 if (debugFn) {
3703 // print polynode and InEle
3704 printf("Ele: %d\n", ele);
3705 printf("ptintsct: %lg, %lg, %lg\n", ptintsct.X, ptintsct.Y,
3706 ptintsct.Z);
3707 printf("nvert: %d\n", nvert);
3708 printf("polynode0: %lg, %lg, %lg\n", polynode[0].X,
3709 polynode[0].Y, polynode[0].Z);
3710 printf("polynode1: %lg, %lg, %lg\n", polynode[1].X,
3711 polynode[1].Y, polynode[1].Z);
3712 printf("polynode2: %lg, %lg, %lg\n", polynode[2].X,
3713 polynode[2].Y, polynode[2].Z);
3714 if (nvert == 4) {
3715 printf("polynode3: %lg, %lg, %lg\n", polynode[3].X,
3716 polynode[3].Y, polynode[3].Z);
3717 }
3718 printf("SumOfAngles: %lg, InEle: %d\n", SumOfAngles, InEle);
3719 fflush(stdout);
3720 }
3721 if (InEle) {
3722 ptintsct.X = (EleArr + ele - 1)->G.Origin.X;
3723 ptintsct.Y = (EleArr + ele - 1)->G.Origin.Y;
3724 ptintsct.Z = (EleArr + ele - 1)->G.Origin.Z;
3725 EleIntsctd = ele;
3726 // Associate this electron to the identified element
3727 NbChUpIonEle[ele]++;
3728 fprintf(fPtIChUpMap, "%d %lg %lg %lg %d %d %d %d\n", ion,
3729 ptintsct.X, ptintsct.Y, ptintsct.Z, prim, InPrim,
3730 ele, InEle);
3731
3732 if (debugFn) {
3733 printf("# ion: %d\n", ion);
3734 printf("%lg %lg %lg\n", xlbend, ylbend, zlbend);
3735 printf("%lg %lg %lg\n", xend, yend, zend);
3736 printf("%lg, %lg, %lg\n", ptintsct.X, ptintsct.Y,
3737 ptintsct.Z);
3738 printf("# Associated primitive: %d\n", prim);
3739 printf(
3740 "# Associated element and origin: %d, %lg, %lg, "
3741 "%lg\n",
3742 ele, (EleArr + ele - 1)->G.Origin.X,
3743 (EleArr + ele - 1)->G.Origin.Y,
3744 (EleArr + ele - 1)->G.Origin.Z);
3745 printf("#NbChUpIonEle on element: %d\n",
3746 NbChUpIonEle[ele]);
3747 fprintf(ftmpIF, "#Element: %d\n", ele);
3748 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[0].X,
3749 polynode[0].Y, polynode[0].Z);
3750 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[1].X,
3751 polynode[1].Y, polynode[1].Z);
3752 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[2].X,
3753 polynode[2].Y, polynode[2].Z);
3754 if (nvert == 4) {
3755 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[3].X,
3756 polynode[3].Y, polynode[3].Z);
3757 }
3758 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[0].X,
3759 polynode[0].Y, polynode[0].Z);
3760 fprintf(ftmpIF, "\n");
3761 fflush(stdout);
3762 }
3763 break; // desired element has been found!
3764 } // if InEle
3765 } // for all elements on this primitive
3766
3767 if (InEle) break;
3769 "Element cannot be identified ... neBEMKnownCharges\n");
3770 return -2;
3771 } // if proper intersection and no extrasection
3772
3773 if ((InPrim) && (intersect) && (!extrasect) && (InEle)) {
3774 // all satisfied
3775 // do not check any further primitive for this electron
3776 break;
3777 }
3778
3779 // If, after checking all the primitives, no interstion is found
3780 // valid
3781 if (prim == (NbPrimitives)) {
3782 // end of the list and no intersection
3783 int nvert;
3784 Point3D polynode[4];
3785 int nearestele = ElementBgn[nearestprim];
3786 double distele = 1.0e6;
3787 double mindistele = 1.0e6; // absurdly high value
3788
3789 if (debugFn) {
3790 printf("prim == (NbPrimitives) ... checking nearest ...\n");
3791 printf("nearestprim: %d, mindist: %lg\n", nearestprim,
3792 mindist);
3793 }
3794
3795 if (mindist <= 10.0e-6) {
3796 PrimIntsctd = nearestprim;
3797 InPrim = 1;
3798 } else {
3799 InPrim = 0;
3800 InEle = 0;
3801 break;
3802 }
3803
3804 for (int ele = ElementBgn[nearestprim]; // check all elements
3805 ele <= ElementEnd[nearestprim]; ++ele) {
3806 nvert = 0;
3807 if ((EleArr + ele - 1)->G.Type == 3) nvert = 3;
3808 if ((EleArr + ele - 1)->G.Type == 4) nvert = 4;
3809 if (!nvert) {
3811 "no vertex element! ... neBEMKnownCharges ...\n");
3812 return -20;
3813 }
3814
3815 /*
3816 polynode[0].X = (EleArr+ele-1)->G.Vertex[0].X;
3817 polynode[0].Y = (EleArr+ele-1)->G.Vertex[0].Y;
3818 polynode[0].Z = (EleArr+ele-1)->G.Vertex[0].Z;
3819 polynode[1].X = (EleArr+ele-1)->G.Vertex[1].X;
3820 polynode[1].Y = (EleArr+ele-1)->G.Vertex[1].Y;
3821 polynode[1].Z = (EleArr+ele-1)->G.Vertex[1].Z;
3822 polynode[2].X = (EleArr+ele-1)->G.Vertex[2].X;
3823 polynode[2].Y = (EleArr+ele-1)->G.Vertex[2].Y;
3824 polynode[2].Z = (EleArr+ele-1)->G.Vertex[2].Z;
3825 if (nvert == 4) {
3826 polynode[3].X = (EleArr+ele-1)->G.Vertex[3].X;
3827 polynode[3].Y = (EleArr+ele-1)->G.Vertex[3].Y;
3828 polynode[3].Z = (EleArr+ele-1)->G.Vertex[3].Z;
3829 }
3830
3831 Vector3D v01, v12, elenorm, unitelenorm;
3832 v01.X = polynode[1].X - polynode[0].X;
3833 v01.Y = polynode[1].Y - polynode[0].Y;
3834 v01.Z = polynode[1].Z - polynode[0].Z;
3835 v12.X = polynode[2].X - polynode[1].X;
3836 v12.Y = polynode[2].Y - polynode[1].Y;
3837 v12.Z = polynode[2].Z - polynode[1].Z;
3838 elenorm = Vector3DCrossProduct(&v01, &v12);
3839 unitelenorm = UnitVector3D(&elenorm);
3840
3841 if ((nvert == 3) || (nvert == 4)) {
3842 if (debugFn) {
3843 printf("nearestprim: %d, element: %d\n",
3844 nearestprim, ele);
3845 printf("vertex0: %lg, %lg, %lg\n", polynode[0].X,
3846 polynode[0].Y, polynode[0].Z);
3847 printf("vertex1: %lg, %lg, %lg\n", polynode[1].X,
3848 polynode[1].Y, polynode[1].Z);
3849 printf("vertex2: %lg, %lg, %lg\n", polynode[2].X,
3850 polynode[2].Y, polynode[2].Z);
3851 if (PrimType[prim] == 4) {
3852 printf("vertex3: %lg, %lg, %lg\n", polynode[3].X,
3853 polynode[3].Y, polynode[3].Z);
3854 }
3855 fprintf(ftmpIF, "#nearestprim: %d, element: %d\n",
3856 nearestprim, ele);
3857 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[0].X,
3858 polynode[0].Y, polynode[0].Z);
3859 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[1].X,
3860 polynode[1].Y, polynode[1].Z);
3861 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[2].X,
3862 polynode[2].Y, polynode[2].Z);
3863 if (PrimType[prim] == 4) {
3864 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[3].X,
3865 polynode[3].Y, polynode[3].Z);
3866 }
3867 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[0].X,
3868 polynode[0].Y, polynode[0].Z);
3869 fprintf(ftmpIF, "\n");
3870 fflush(stdout);
3871 } // debugFn
3872
3873 // use a, b, c (normal is ai + bj + ck)
3874 // at one of the nodes to get d
3875 // ax + by + cz + d = 0 is the equation of the plane
3876 double a = unitelenorm.X;
3877 double b = unitelenorm.Y;
3878 double c = unitelenorm.Z;
3879 double d = - a * polynode[0].X - b * polynode[0].Y - c * polynode[0].Z;
3880 // distance of the end point to this primitve is
3881 distele = (xend * a + yend * b + zend * c + d) /
3882 sqrt(a * a + b * b + c * c);
3883 distele = fabs(distele); // if only magnitude is required
3884 */
3885
3886 Vector3D eleOrigin;
3887 eleOrigin.X = (EleArr + ele - 1)->G.Origin.X;
3888 eleOrigin.Y = (EleArr + ele - 1)->G.Origin.Y;
3889 eleOrigin.Z = (EleArr + ele - 1)->G.Origin.Z;
3890 distele = (eleOrigin.X - xend) * (eleOrigin.X - xend) +
3891 (eleOrigin.Y - yend) * (eleOrigin.Y - yend) +
3892 (eleOrigin.Z - zend) * (eleOrigin.Z - zend);
3893 distele = pow(distele, 0.5);
3894
3895 if (ele == ElementBgn[nearestprim]) {
3896 mindistele = distele;
3897 nearestele = ele;
3898 }
3899 if (distele < mindistele) {
3900 mindistele = distele;
3901 nearestele = ele;
3902 }
3903
3904 if (debugFn) {
3905 // printf("a, b, c, d, dist: %lg, %lg, %lg, %lg, %lg\n",
3906 // a, b, c, d, dist);
3907 // printf("vector n: ai + bj + ck\n");
3908 // printf("vector v: xgrd, ygrd, zgrd: %lg, %lg, %lg\n",
3909 // xgrd, ygrd, zgrd);
3910 printf(
3911 "distele: %lg, mindist: %lg, from nearest ele: %d\n",
3912 distele, mindistele, nearestele);
3913 fflush(stdout);
3914 }
3915
3916 // } // if PrimType is 3 or 4
3917 } // for elements in nearestprim
3918
3919 // if(mindistele <= 10.0e-6)
3920 // {
3921 EleIntsctd = nearestele;
3922 InEle = 1;
3923 ptintsct.X = (EleArr + EleIntsctd - 1)->G.Origin.X;
3924 ptintsct.Y = (EleArr + EleIntsctd - 1)->G.Origin.Y;
3925 ptintsct.Z = (EleArr + EleIntsctd - 1)->G.Origin.Z;
3926 NbChUpIonEle[EleIntsctd]++;
3927
3928 fprintf(fPtIChUpMap, "%d %lg %lg %lg %d %d %d %d\n", ion,
3929 ptintsct.X, ptintsct.Y, ptintsct.Z, PrimIntsctd, InPrim,
3930 EleIntsctd, InEle);
3931 // }
3932
3933 if (debugFn) {
3934 printf("# ion: %d\n", ion);
3935 printf("%lg %lg %lg\n", xlbend, ylbend, zlbend);
3936 printf("%lg %lg %lg\n", xend, yend, zend);
3937 printf("%lg, %lg, %lg\n", ptintsct.X, ptintsct.Y, ptintsct.Z);
3938 printf("# Associated primitive: %d\n", PrimIntsctd);
3939 printf("# Associated element and origin: %d, %lg, %lg, %lg\n",
3940 EleIntsctd, (EleArr + EleIntsctd - 1)->G.Origin.X,
3941 (EleArr + EleIntsctd - 1)->G.Origin.Y,
3942 (EleArr + EleIntsctd - 1)->G.Origin.Z);
3943 printf("#NbChUpIonEle on element: %d\n",
3944 NbChUpIonEle[EleIntsctd]);
3945 fprintf(ftmpIF, "#Element: %d\n", EleIntsctd);
3946 polynode[0].X = (EleArr + EleIntsctd - 1)->G.Vertex[0].X;
3947 polynode[0].Y = (EleArr + EleIntsctd - 1)->G.Vertex[0].Y;
3948 polynode[0].Z = (EleArr + EleIntsctd - 1)->G.Vertex[0].Z;
3949 polynode[1].X = (EleArr + EleIntsctd - 1)->G.Vertex[1].X;
3950 polynode[1].Y = (EleArr + EleIntsctd - 1)->G.Vertex[1].Y;
3951 polynode[1].Z = (EleArr + EleIntsctd - 1)->G.Vertex[1].Z;
3952 polynode[2].X = (EleArr + EleIntsctd - 1)->G.Vertex[2].X;
3953 polynode[2].Y = (EleArr + EleIntsctd - 1)->G.Vertex[2].Y;
3954 polynode[2].Z = (EleArr + EleIntsctd - 1)->G.Vertex[2].Z;
3955 if (nvert == 4) {
3956 polynode[3].X = (EleArr + EleIntsctd - 1)->G.Vertex[3].X;
3957 polynode[3].Y = (EleArr + EleIntsctd - 1)->G.Vertex[3].Y;
3958 polynode[3].Z = (EleArr + EleIntsctd - 1)->G.Vertex[3].Z;
3959 }
3960 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[0].X, polynode[0].Y,
3961 polynode[0].Z);
3962 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[1].X, polynode[1].Y,
3963 polynode[1].Z);
3964 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[2].X, polynode[2].Y,
3965 polynode[2].Z);
3966 if (nvert == 4) {
3967 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[3].X,
3968 polynode[3].Y, polynode[3].Z);
3969 }
3970 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[0].X, polynode[0].Y,
3971 polynode[0].Z);
3972 fprintf(ftmpIF, "\n");
3973 fflush(stdout);
3974 } // debug
3975 } // if prim == NbPrimitives
3976
3977 } // for all primitives // just not those on the volume
3978
3979 if (debugFn) // check ion positions, volume primitives and elements
3980 {
3981 char ionposdbg[256], inbstr[10];
3982 sprintf(inbstr, "%d", ion);
3983 strcpy(ionposdbg, "/tmp/Ion");
3984 strcat(ionposdbg, inbstr);
3985 strcat(ionposdbg, ".out");
3986 FILE *fipd = fopen(ionposdbg, "w");
3987 if (fipd == NULL) {
3988 printf(
3989 "cannot open writable file to debug ion positions ...\n");
3990 printf("returning ...\n");
3991 return -111;
3992 }
3993 // write electron number, end, lbend, volume, primitive, elements,
3994 // intxn
3995 fprintf(fipd, "#ion: %d %d\n", inb, ion); // should print same
3996 fprintf(fipd, "#last but end position:\n");
3997 fprintf(fipd, "%lg %lg %lg\n", xlbend, ylbend, zlbend);
3998 fprintf(fipd, "#end position:\n");
3999 fprintf(fipd, "%lg %lg %lg\n\n", xend, yend, zend);
4000
4001 fprintf(fipd, "#intersected primitive number: %d\n", PrimIntsctd);
4002 if (PrimIntsctd >= 1) {
4003 fprintf(fipd, "#PrimType: %d\n", PrimType[PrimIntsctd]);
4004 fprintf(fipd, "#prim vertices:\n");
4005 fprintf(fipd, "%lg %lg %lg\n", XVertex[PrimIntsctd][0],
4006 YVertex[PrimIntsctd][0], ZVertex[PrimIntsctd][0]);
4007 fprintf(fipd, "%lg %lg %lg\n", XVertex[PrimIntsctd][1],
4008 YVertex[PrimIntsctd][1], ZVertex[PrimIntsctd][1]);
4009 fprintf(fipd, "%lg %lg %lg\n", XVertex[PrimIntsctd][2],
4010 YVertex[PrimIntsctd][2], ZVertex[PrimIntsctd][2]);
4011 if (PrimType[PrimIntsctd] == 4) {
4012 fprintf(fipd, "%lg %lg %lg\n", XVertex[PrimIntsctd][3],
4013 YVertex[PrimIntsctd][3], ZVertex[PrimIntsctd][3]);
4014 }
4015 fprintf(fipd, "%lg %lg %lg\n", XVertex[PrimIntsctd][0],
4016 YVertex[PrimIntsctd][0], ZVertex[PrimIntsctd][0]);
4017 fprintf(fipd, "\n");
4018
4019 fprintf(fipd, "#ptintsct:\n");
4020 fprintf(fipd, "%lg %lg %lg\n", ptintsct.X, ptintsct.Y,
4021 ptintsct.Z);
4022 fprintf(fipd, "\n");
4023 }
4024
4025 fprintf(fipd, "#intersected element number: %d\n", EleIntsctd);
4026 if (EleIntsctd >= 1) {
4027 int gtype = (EleArr + EleIntsctd - 1)->G.Type;
4028 fprintf(fipd, "#EleType: %d\n", gtype);
4029 fprintf(fipd, "#element vertices:\n");
4030 double x0 = (EleArr + EleIntsctd - 1)->G.Vertex[0].X;
4031 double y0 = (EleArr + EleIntsctd - 1)->G.Vertex[0].Y;
4032 double z0 = (EleArr + EleIntsctd - 1)->G.Vertex[0].Z;
4033 double x1 = (EleArr + EleIntsctd - 1)->G.Vertex[1].X;
4034 double y1 = (EleArr + EleIntsctd - 1)->G.Vertex[1].Y;
4035 double z1 = (EleArr + EleIntsctd - 1)->G.Vertex[1].Z;
4036 double x2 = (EleArr + EleIntsctd - 1)->G.Vertex[2].X;
4037 double y2 = (EleArr + EleIntsctd - 1)->G.Vertex[2].Y;
4038 double z2 = (EleArr + EleIntsctd - 1)->G.Vertex[2].Z;
4039 fprintf(fipd, "%lg %lg %lg\n", x0, y0, z0);
4040 fprintf(fipd, "%lg %lg %lg\n", x1, y1, z1);
4041 fprintf(fipd, "%lg %lg %lg\n", x2, y2, z2);
4042 if (gtype == 4) {
4043 double x3 = (EleArr + EleIntsctd - 1)->G.Vertex[3].X;
4044 double y3 = (EleArr + EleIntsctd - 1)->G.Vertex[3].Y;
4045 double z3 = (EleArr + EleIntsctd - 1)->G.Vertex[3].Z;
4046 fprintf(fipd, "%lg %lg %lg\n", x3, y3, z3);
4047 }
4048 fprintf(fipd, "%lg %lg %lg\n", x0, y0, z0);
4049 fprintf(fipd, "\n");
4050
4051 fprintf(fipd, "#ptintsct:\n");
4052 fprintf(fipd, "%lg %lg %lg\n", ptintsct.X, ptintsct.Y,
4053 ptintsct.Z);
4054 fprintf(fipd, "\n");
4055 }
4056 fclose(fipd);
4057 } // if 1
4058 } // for all the ions
4059 fclose(fPtIChUpMap);
4060
4061 // This file contains information about number of ions (I)
4062 // and total (E+I) charge deposition on each element
4063 FILE *fEleEIChUpMap = fopen("EleE+IChUpMap.out", "w");
4064 if (fEleEIChUpMap == NULL) {
4065 printf("cannot open EleE+IChUpMap.out file for writing ...\n");
4066 return 111;
4067 }
4068 for (int ele = 1; ele <= NbElements; ++ele) {
4069 (EleArr + ele - 1)->Assigned +=
4070 ChUpFactor * Q_I * NbChUpIonEle[ele] / (EleArr + ele - 1)->G.dA;
4071 fprintf(fEleEIChUpMap, "%d %lg %lg %lg %d %lg\n", ele,
4072 (EleArr + ele - 1)->G.Origin.X,
4073 (EleArr + ele - 1)->G.Origin.Y,
4074 (EleArr + ele - 1)->G.Origin.Z, NbChUpIonEle[ele],
4075 (EleArr + ele - 1)->Assigned);
4076 }
4077 fclose(fEleEIChUpMap);
4078
4079 fclose(ftmpIF);
4080 free(NbChUpIonEle);
4081 } // Calculation for ions ends
4082
4083 } // OptChargingUp
4084
4085 fclose(ChargingUpInpFile);
4086 } // else ChargingUpInpFile
4087
4088 if (debugFn) {
4089 // print all the elements and their number of charging up e-s and i-s
4090 }
4091 } // charging up parameters set up
4092
4093 return (0);
4094} // neBEMChargingUp ends
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:337
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314
int neBEMGetNbOfLines(const char fname[])
double neBEMChkInPoly(int n, Point3D *p, Point3D q)
#define neBEMtwopi
neBEMGLOBAL Element * EleArr
Definition: neBEM.h:169
neBEMGLOBAL double * ZNorm
Definition: neBEM.h:71
#define Q_E
Definition: neBEM.h:27
neBEMGLOBAL double ** ZVertex
Definition: neBEM.h:71
neBEMGLOBAL int OptChargingUp
Definition: neBEM.h:53
neBEMGLOBAL double * XNorm
Definition: neBEM.h:71
neBEMGLOBAL int NbPrimitives
Definition: neBEM.h:57
neBEMGLOBAL int NbElements
Definition: neBEM.h:113
neBEMGLOBAL double ** YVertex
Definition: neBEM.h:71
#define Q_I
Definition: neBEM.h:28
neBEMGLOBAL double * YNorm
Definition: neBEM.h:71
neBEMGLOBAL double ** XVertex
Definition: neBEM.h:71
neBEMGLOBAL int * ElementEnd
Definition: neBEM.h:94
neBEMGLOBAL int * InterfaceType
Definition: neBEM.h:68
neBEMGLOBAL int * PrimType
Definition: neBEM.h:67
neBEMGLOBAL int * ElementBgn
Definition: neBEM.h:94
Definition: Vector.h:21
double X
Definition: Vector.h:22
double Z
Definition: Vector.h:24
double Y
Definition: Vector.h:23
double Z
Definition: Vector.h:31
double Y
Definition: Vector.h:30
double X
Definition: Vector.h:29

Referenced by ComputeSolution().

◆ neBEMChkInPoly()

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

◆ neBEMDeleteAllWeightingFields()

INTFACEGLOBAL void neBEMDeleteAllWeightingFields ( void  )

Definition at line 4666 of file neBEMInterface.c.

4666 {
4667 const int MaxWtField = 100; // being used while allocating memory
4668 for (int id = 1; id < MaxWtField; ++id) { // count from 1
4669 free(WtFieldChDen[id]);
4670 free(AvWtChDen[id]);
4671 }
4672 free(WtFieldChDen);
4673 free(AvWtChDen);
4674}
neBEMGLOBAL double ** AvWtChDen
Definition: neBEM.h:333
neBEMGLOBAL double ** WtFieldChDen
Definition: neBEM.h:333

◆ neBEMDeleteWeightingField()

INTFACEGLOBAL void neBEMDeleteWeightingField ( int  IdWtField)

Definition at line 4660 of file neBEMInterface.c.

4660 {
4661 free(WtFieldChDen[IdWtField]);
4662 free(AvWtChDen[IdWtField]);
4663}

◆ neBEMDiscretize()

INTFACEGLOBAL int neBEMDiscretize ( int **  elementNbs)

Definition at line 1913 of file neBEMInterface.c.

1913 {
1914 // int fstatus;
1915
1916 // For a model and a mesh that were defined before and for which data were
1917 // stored in two files - one for primitives and one for elements
1918 // The following operation does not assume that the primitives have been read
1919 // in from a stored file. In essence, these two read-in-s are maintained
1920 // independent of each other. However, for a stored element file to be useful,
1921 // both the model and mesh have to be old (not new).
1922 if ((!NewModel) && (!NewMesh) && (!NewBC) && (OptStoreElements)) {
1923 int fstatus = ReadElements();
1924 if (fstatus) {
1925 neBEMMessage("neBEMDiscretize - problem reading stored Elements.\n");
1926 return -1;
1927 }
1928 neBEMState = 4;
1929 return 0;
1930 }
1931
1932 // Otherwise, continue with fresh discretization
1933 if (neBEMState != 3) {
1934 printf("discretization can continue only in State 3 ...\n");
1935 return -1;
1936 }
1937
1938 // Only the number of primitives has been ascertained.
1939 // All the rest globally important numbers will be determined in this function
1940 NbSurfs = 0;
1941 NbWires = 0;
1942 NbElements = 0;
1943
1944 // Here, the primitive can be analyzed and the elements necessary to
1945 // discretize it, may be determined. A user hint may help that can be supplied
1946 // during the setting up of the device. The hint may be as naive as gross,
1947 // coarse, medium, regular, fine depending on which the element size (in %
1948 // of the device / primitive) may be decided upon.
1949 // Earlier method of specifying the number of primitives on a priori has been
1950 // over-ridden.
1951 // Now we prescribe the length / area of each element and discretize each
1952 // primitive such that each element has an length / area close to the
1953 // suggested value.
1954 // Since the number of elements are being decided here, all the shapes will
1955 // have to be one among wire, right triangle and restangle. Any decomposition
1956 // of arbitrary polygons into rectangles (as many as possible) and right
1957 // triangles will have to be done earlier. All the counts have been
1958 // incremented by one to be on the safe side! The additional memory can be
1959 // minimized through a more careful computation of the required allocation for
1960 // each type of primitive.
1961 char MeshLogFile[256];
1962
1963 strcpy(MeshLogFile, MeshOutDir);
1964 strcat(MeshLogFile, "/MeshLog.out");
1965 fMeshLog = fopen(MeshLogFile, "w");
1966 fprintf(fMeshLog, "Details of primitive discretization\n");
1967
1968 for (int prim = 1; prim <= NbPrimitives; ++prim) {
1969 if (NbVertices[prim] == 4) {
1970 NbSurfSegX[prim] = NbElemsOnPrimitives[prim][1];
1971 NbSurfSegZ[prim] = NbElemsOnPrimitives[prim][2];
1972 int fstatus = AnalyzePrimitive(prim, &NbSurfSegX[prim], &NbSurfSegZ[prim]);
1973 if (fstatus == 0) {
1974 neBEMMessage("neBEMDiscretize - AnalyzePrimitve");
1975 return -1;
1976 }
1977 NbElements += (NbSurfSegX[prim] + 1) * (NbSurfSegZ[prim] + 1);
1978 }
1979 if (NbVertices[prim] == 3) {
1980 NbSurfSegX[prim] = NbElemsOnPrimitives[prim][1];
1981 NbSurfSegZ[prim] = NbElemsOnPrimitives[prim][2];
1982 int fstatus = AnalyzePrimitive(prim, &NbSurfSegX[prim], &NbSurfSegZ[prim]);
1983 if (fstatus == 0) {
1984 neBEMMessage("neBEMDiscretize - AnalyzePrimitive");
1985 return -1;
1986 }
1987 NbElements += (NbSurfSegX[prim] + 1) * (NbSurfSegZ[prim] + 1);
1988 }
1989 if (NbVertices[prim] == 2) {
1990 int itmp;
1991 NbWireSeg[prim] = NbElemsOnPrimitives[prim][1];
1992 int fstatus = AnalyzePrimitive(prim, &NbWireSeg[prim], &itmp);
1993 if (fstatus == 0) {
1994 neBEMMessage("neBEMDiscretize - AnalyzePrimitive");
1995 return -1;
1996 }
1997 NbElements += (NbWireSeg[prim] + 1);
1998 }
1999 if (DebugLevel == 101) {
2000 if (NbVertices[prim] == 2) {
2001 printf("Primitive %d to be discretized into %d elements.\n", prim,
2002 NbWireSeg[prim]);
2003 } else {
2004 printf("Primitive %d to be discretized into %d X %d elements.\n", prim,
2005 NbSurfSegX[prim], NbSurfSegZ[prim]);
2006 }
2007 }
2008 }
2009 printf("Memory allocated for maximum %d elements.\n", NbElements);
2010 fclose(fMeshLog);
2011
2012 // Allocate enough space to store all the elements
2013 if (neBEMState == 3) {
2014 printf("neBEMDiscretize: NbElements = %d, sizeof(Element) = %zu\n",
2015 NbElements, sizeof(Element));
2016 if (EleArr) {
2017 Element *tmp = (Element *)realloc(EleArr, NbElements * sizeof(Element));
2018 if (tmp != NULL) {
2019 EleArr = tmp;
2020 EleCntr = 0;
2021 } else {
2022 free(EleArr);
2023 printf("neBEMDiscretize: Re-allocating EleArr failed.\n");
2024 return (1);
2025 }
2026 printf("neBEMDiscretize: Re-allocated EleArr.\n");
2027 } // if EleArr => re-allocation
2028 else {
2029 EleArr = (Element *)malloc(NbElements * sizeof(Element));
2030 if (EleArr == NULL) {
2031 neBEMMessage("neBEMDiscretize - EleArr malloc");
2032 return -1;
2033 }
2034 } // else EleArr => fresh allocation
2035 } // neBEMState == 3
2036
2037 // Prepare a data file that will contain the plotting information of the
2038 // the primitives and the elements
2039 if (OptGnuplot) {
2040 char GnuFile[256];
2041 strcpy(GnuFile, MeshOutDir);
2042 strcat(GnuFile, "/GViewDir/gPrimView.gp");
2043 fgnuPrim = fopen(GnuFile, "w");
2044 fprintf(fgnuPrim, "set title \"neBEM primitives in gnuplot VIEWER\"\n");
2045 // fprintf(fgnu, "#set label 1 \'LengthScale = %d\', LengthScale, right\n");
2046 fprintf(fgnuPrim, "#set pm3d\n");
2047 fprintf(fgnuPrim, "#set style data pm3d\n");
2048 fprintf(fgnuPrim, "#set palette model CMY\n");
2049 fprintf(fgnuPrim, "set hidden3d\n");
2050 fprintf(fgnuPrim, "set nokey\n");
2051 fprintf(fgnuPrim, "set xlabel \"X\"\n");
2052 fprintf(fgnuPrim, "set ylabel \"Y\"\n");
2053 fprintf(fgnuPrim, "set zlabel \"Z\"\n");
2054 fprintf(fgnuPrim, "set view 70, 335, 1, 1\n");
2055 fprintf(fgnuPrim, "\nsplot \\\n");
2056
2057 strcpy(GnuFile, MeshOutDir);
2058 strcat(GnuFile, "/GViewDir/gElemView.gp");
2059 fgnuElem = fopen(GnuFile, "w");
2060 fprintf(fgnuElem, "set title \"neBEM elements in gnuplot VIEWER\"\n");
2061 // fprintf(fgnu, "#set label 1 \'LengthScale = %d\', LengthScale, right\n");
2062 fprintf(fgnuElem, "#set pm3d\n");
2063 fprintf(fgnuElem, "#set style data pm3d\n");
2064 fprintf(fgnuElem, "#set palette model CMY\n");
2065 fprintf(fgnuElem, "set hidden3d\n");
2066 fprintf(fgnuElem, "set nokey\n");
2067 fprintf(fgnuElem, "set xlabel \"X\"\n");
2068 fprintf(fgnuElem, "set ylabel \"Y\"\n");
2069 fprintf(fgnuElem, "set zlabel \"Z\"\n");
2070 fprintf(fgnuElem, "set view 70, 335, 1, 1\n");
2071 fprintf(fgnuElem, "\nsplot \\\n");
2072
2073 strcpy(GnuFile, MeshOutDir);
2074 strcat(GnuFile, "/GViewDir/gMeshView.gp");
2075 fgnuMesh = fopen(GnuFile, "w");
2076 fprintf(fgnuMesh, "set title \"neBEM mesh in gnuplot VIEWER\"\n");
2077 // fprintf(fgnu, "#set label 1 \'LengthScale = %d\', LengthScale, right\n");
2078 fprintf(fgnuMesh, "#set pm3d\n");
2079 fprintf(fgnuMesh, "#set style data pm3d\n");
2080 fprintf(fgnuMesh, "#set palette model CMY\n");
2081 fprintf(fgnuMesh, "set hidden3d\n");
2082 fprintf(fgnuMesh, "set nokey\n");
2083 fprintf(fgnuMesh, "set xlabel \"X\"\n");
2084 fprintf(fgnuMesh, "set ylabel \"Y\"\n");
2085 fprintf(fgnuMesh, "set zlabel \"Z\"\n");
2086 fprintf(fgnuMesh, "set view 70, 335, 1, 1\n");
2087 fprintf(fgnuMesh, "\nsplot \\\n");
2088 }
2089
2090 for (int prim = 1; prim <= NbPrimitives; ++prim) {
2091 switch (PrimType[prim]) {
2092 int fstatus;
2093 case 3: // triangular surface
2094 case 4: // rectangular surface
2095 ++NbSurfs;
2096 fstatus = SurfaceElements(
2097 prim, NbVertices[prim], XVertex[prim], YVertex[prim], ZVertex[prim],
2098 XNorm[prim], YNorm[prim], ZNorm[prim], VolRef1[prim], VolRef2[prim],
2099 InterfaceType[prim], ApplPot[prim], ApplCh[prim], Lambda[prim],
2100 NbSurfSegX[prim], NbSurfSegZ[prim]);
2101 if (fstatus != 0) {
2102 neBEMMessage("neBEMDiscretize - SurfaceElements");
2103 return -1;
2104 }
2105 break;
2106 case 2: // wire - a wire presumably has only 2 vertices
2107 ++NbWires; // it has one radius and one segmentation information
2108 fstatus = WireElements(
2109 prim, NbVertices[prim], XVertex[prim], YVertex[prim], ZVertex[prim],
2110 Radius[prim], VolRef1[prim], VolRef2[prim], InterfaceType[prim],
2111 ApplPot[prim], ApplCh[prim], Lambda[prim], NbWireSeg[prim]);
2112 if (fstatus != 0) {
2113 neBEMMessage("neBEMDiscretize - WireElements");
2114 return -1;
2115 }
2116 break;
2117
2118 default:
2119 printf("PrimType out of range in CreateElements ... exiting ...\n");
2120 exit(-1);
2121 } // switch PrimType ends
2122 } // loop on prim number ends
2123
2124 if (OptGnuplot) {
2125 fprintf(fgnuPrim, "\n\npause-1");
2126 fclose(fgnuPrim);
2127 fprintf(fgnuElem, "\n\npause-1");
2128 fclose(fgnuElem);
2129 fprintf(fgnuMesh, "\n\npause-1");
2130 fclose(fgnuMesh);
2131 }
2132
2133 // If the required memory exceeds the maximum allowed number of elements
2134 if (EleCntr > NbElements) {
2135 neBEMMessage("neBEMDiscretize - EleCntr more than NbElements!");
2136 return -1;
2137 }
2138
2139 // Check whether collocation points overlap
2140 {
2141 int startcntr = 1, cntr1, cntr2, NbCollPtOverlaps = 0;
2142 Point3D pt1, pt2;
2143 double dist;
2144 for (cntr1 = startcntr; cntr1 <= EleCntr; ++cntr1) {
2145 pt1.X = (EleArr + cntr1 - 1)->BC.CollPt.X;
2146 pt1.Y = (EleArr + cntr1 - 1)->BC.CollPt.Y;
2147 pt1.Z = (EleArr + cntr1 - 1)->BC.CollPt.Z;
2148 for (cntr2 = cntr1 + 1; cntr2 <= EleCntr; ++cntr2) {
2149 pt2.X = (EleArr + cntr2 - 1)->BC.CollPt.X;
2150 pt2.Y = (EleArr + cntr2 - 1)->BC.CollPt.Y;
2151 pt2.Z = (EleArr + cntr2 - 1)->BC.CollPt.Z;
2152
2153 dist = GetDistancePoint3D(&pt1, &pt2);
2154 if (dist <=
2155 MINDIST) // we need a linked-list here so that the overlapped
2156 { // element is easily deleted and the rest upgraded immediately
2157 ++NbCollPtOverlaps;
2158
2159 // Upgrade the element array manually, starting from cntr2 and restart
2160 // the overlap check. At present it is only a warning to the user with
2161 // some relevant information.
2162 // Find the primitives and volumes for the overlapping elements
2163 // The element structure should also maintain information on the
2164 // volumes that an element belongs to.
2165 int prim1 = (EleArr + cntr1 - 1)->PrimitiveNb;
2166 int volele1 = VolRef1[prim1];
2167 int prim2 = (EleArr + cntr2 - 1)->PrimitiveNb;
2168 int volele2 = VolRef1[prim2];
2169
2170 neBEMMessage("neBEMDiscretize - Overlapping collocation points!");
2171 printf("Element %d, primitive %d, volume %d overlaps with\n", cntr1,
2172 prim1, volele1);
2173 printf("\telement %d, primitive %d, volume %d.\n", cntr2, prim2,
2174 volele2);
2175 printf("\tposition 1: (%g , %g , %g) micron,\n", 1e6 * pt1.X,
2176 1e6 * pt1.Y, 1e6 * pt1.Z);
2177 printf("\tposition 2: (%g , %g , %g) micron.\n", 1e6 * pt2.X,
2178 1e6 * pt2.Y, 1e6 * pt2.Z);
2179 printf("Please redo the geometry.\n");
2180 return -1;
2181 } // if dist <= MINDIST
2182 } // for cntr2
2183 } // for cntr1
2184 } // check collocation point overlap
2185
2186 NbElements = EleCntr; // the final number of elements
2187 printf("Total final number of elements: %d\n", NbElements);
2188
2189 // Store element related data in a file for a new mesh created, if opted for
2190 if (NewMesh && OptStoreElements) {
2191 if (OptFormattedFile) {
2192 int fstatus = WriteElements();
2193 if (fstatus) {
2194 neBEMMessage("neBEMDiscretize - problem writing Elements.\n");
2195 return -1;
2196 }
2197 } // formatted file
2198
2199 if (OptUnformattedFile) {
2201 "neBEMDiscretize - unformatted write not inplemented yet.\n");
2202 return -1;
2203 } // unformatted file
2204 } // store elements
2205
2206 neBEMState = 4;
2207 stopClock = clock();
2209 printf("to complete discretization\n");
2210
2211 return (0);
2212} // neBEMDiscretize ends
#define MINDIST
Definition: Isles.h:17
double GetDistancePoint3D(Point3D *a, Point3D *b)
Definition: Vector.c:192
int ReadElements(void)
int WriteElements(void)
INTFACEGLOBAL FILE * fgnuPrim
INTFACEGLOBAL FILE * fgnuElem
INTFACEGLOBAL FILE * fgnuMesh
INTFACEGLOBAL int OptGnuplot
neBEMGLOBAL double * ApplPot
Definition: neBEM.h:75
neBEMGLOBAL int DebugLevel
Definition: neBEM.h:235
neBEMGLOBAL int EleCntr
Definition: neBEM.h:112
neBEMGLOBAL int * NbWireSeg
Definition: neBEM.h:99
neBEMGLOBAL int * NbSurfSegX
Definition: neBEM.h:98
neBEMGLOBAL int OptFormattedFile
Definition: neBEM.h:49
neBEMGLOBAL int * NbVertices
Definition: neBEM.h:70
neBEMGLOBAL int WireElements(int prim, int nvertex, double xvert[], double yvert[], double zvert[], double radius, int volref1, int volref2, int inttype, double potential, double charge, double lambda, int NbWireSeg)
Definition: ReTriM.c:82
neBEMGLOBAL int AnalyzePrimitive(int, int *, int *)
Definition: ReTriM.c:121
neBEMGLOBAL double * Lambda
Definition: neBEM.h:75
neBEMGLOBAL int * NbSurfSegZ
Definition: neBEM.h:98
neBEMGLOBAL int NbWires
Definition: neBEM.h:97
neBEMGLOBAL int * VolRef1
Definition: neBEM.h:76
neBEMGLOBAL int * VolRef2
Definition: neBEM.h:76
neBEMGLOBAL int NbSurfs
Definition: neBEM.h:97
neBEMGLOBAL double * Radius
Definition: neBEM.h:72
neBEMGLOBAL int OptUnformattedFile
Definition: neBEM.h:50
neBEMGLOBAL FILE * fMeshLog
Definition: neBEM.h:114
neBEMGLOBAL int SurfaceElements(int prim, int nvertex, double xvert[], double yvert[], double zvert[], double xnorm, double ynorm, double znorm, int volref1, int volref2, int inttype, double potential, double charge, double lambda, int NbSegX, int NbSegZ)
Definition: ReTriM.c:33
neBEMGLOBAL int OptStoreElements
Definition: neBEM.h:46
neBEMGLOBAL double * ApplCh
Definition: neBEM.h:75
Definition: neBEM.h:155

◆ neBEMEnd()

INTFACEGLOBAL int neBEMEnd ( void  )

Definition at line 4756 of file neBEMInterface.c.

4756 {
4757 fprintf(fIsles,
4758 "IslesCntr: %d, ExactCntr: %d, FailureCntr: %d, ApproxCntr: %d\n",
4760 fclose(fIsles);
4761 fIsles = NULL;
4762 printf("neBEM ends ... bye!\n");
4763
4764 return 0;
4765} // neBEMEnd ends
ISLESGLOBAL int ApproxCntr
Definition: Isles.h:34
ISLESGLOBAL int IslesCntr
Definition: Isles.h:34
ISLESGLOBAL FILE * fIsles
Definition: Isles.h:36
ISLESGLOBAL int FailureCntr
Definition: Isles.h:34
ISLESGLOBAL int ExactCntr
Definition: Isles.h:34

◆ neBEMField()

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

Definition at line 4561 of file neBEMInterface.c.

4561 {
4562 if (neBEMState < 9) {
4563 printf("neBEMField cannot be called before reaching state 9.\n");
4564 return (-1);
4565 }
4566
4567 // printf("neBEMField called %8d times", ++neBEMFieldCallCntr);
4568
4569 double Pot;
4570 int fstatus;
4571 if (OptFastVol) // Note: this is not the Create or Read option
4572 {
4573 fstatus = FastPFAtPoint(point, &Pot, field);
4574 if (fstatus != 0) {
4575 neBEMMessage("neBEMField - FastPFAtPoint");
4576 return -1;
4577 }
4578 } else {
4579 fstatus = PFAtPoint(point, &Pot, field);
4580 if (fstatus != 0) {
4581 neBEMMessage("neBEMField - PFAtPoint");
4582 return -1;
4583 }
4584 }
4585
4586 *potential = Pot;
4587
4588 // 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");
4589
4590 return (0);
4591} // neBEMField ends
int FastPFAtPoint(Point3D *globalP, double *Potential, Vector3D *globalF)
int PFAtPoint(Point3D *globalP, double *Potential, Vector3D *globalF)
neBEMGLOBAL int OptFastVol
Definition: neBEM.h:403

◆ 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.

39 {
40
41 // Version information
42 strcpy(neBEMVersion, "1.9.08");
43 strcpy(ISLESVersion, "1.4.8");
44 printf("Using neBEM version %s and ISLES version %s\n", neBEMVersion,
46
47 // The following is not an absolute necessity, but can ease the burden of the
48 // user by setting default values of some of the important variable. Please
49 // take a look at the src/Interface/ExampleDev2neBEM.c to find out what values
50 // are being set, and what they mean.
51 int fstatus = neBEMSetDefaults();
52 if (fstatus != 0) {
53 neBEMMessage("neBEMInitialize - neBEMSetDefaults");
54 return -1;
55 }
56
57 // Change some of the global variables according to requirements.
58 // Note that the solution flags and counters are set before the neBEMSolve
59 // is invoked.
60 LengthScale = 1.0;
61 DebugLevel = 0;
62
63 // The following function allows input files to be read and the geometry
64 // set according to these files, as was customary in the previous versions
65 if (OptDeviceFile) {
66 printf("Reading geometry details from %s\n", DeviceInputFile);
67 fstatus = neBEMGetInputsFromFiles();
68 if (fstatus != 0) {
69 neBEMMessage("neBEMInitialize - neBEMGetInputFromFiles");
70 return -1;
71 }
72 }
73
74 // creation of the mother output directory should be necessary only once
75 if (neBEMState == 0) {
76 fstatus = CreateDirStr();
77 if (fstatus != 0) {
78 neBEMMessage("neBEMInitialize - CreateDirStr");
79 return -1;
80 }
81 }
82
83 // Create Isles log file for keeping track of approximations (numerical
84 // quadrature) in case evaluation of algebraic expressions fails.
85 char IslesFile[256];
86 strcpy(IslesFile, PPOutDir);
87 strcat(IslesFile, "/Isles.log");
88 fIsles = fopen(IslesFile, "w");
89 if (fIsles == NULL) {
90 neBEMMessage("neBEMInitialize - IslesFile");
91 return -1;
92 }
93
94 // following integers keep track of the success / failure of the exact
95 // expressions in estimating realistic physical properties.
97
98 // Set up parameters related to neBEM computations
99 int RqstdThreads = 1;
100 FILE *processFile = fopen("neBEMProcess.inp", "r");
101 if (processFile == NULL) {
102 printf("neBEMProcess.inp absent ... assuming defaults ...\n");
103 PrimAfter = 0;
104 if (NbThreads > 0) RqstdThreads = NbThreads;
105 } else {
106 fscanf(processFile, "PrimAfter: %d\n", &PrimAfter);
107 fscanf(processFile, "RqstdThreads: %d\n", &RqstdThreads);
108 fclose(processFile);
109 }
110
111#ifdef _OPENMP
112 int MaxProcessors = omp_get_num_procs();
113 if (RqstdThreads > 1) {
114 if (RqstdThreads < MaxProcessors) {
115 // one processor left alone
116 omp_set_num_threads(RqstdThreads);
117 } else {
118 printf("RqstdThreads: %d\n", RqstdThreads);
119 RqstdThreads = MaxProcessors - 1;
120 omp_set_num_threads(RqstdThreads);
121 printf("Adjusted RqstdThreads: %d\n", RqstdThreads);
122 }
123 } else {
124 // Work with one thread
125 RqstdThreads = 1; // cannot be zero or negative!
126 omp_set_num_threads(RqstdThreads);
127 printf("RqstdThreads: %d => No Multi-threading ...\n", RqstdThreads);
128 }
129
130 // OpenMP related information
131 printf("PrimAfter: %d\n", PrimAfter);
132 printf("RqstdThreads: %d, MaxProcessors: %d\n", RqstdThreads, MaxProcessors);
133 printf("Maximum number of threads to be used for parallelization: %d\n",
134 omp_get_max_threads());
135 printf("Number of threads used for neBEMInitialize: %d\n", omp_get_num_threads());
136#endif
137 // Set up parameters related to voxelized data export for Garfield++
138 FILE *voxelInpFile = fopen("neBEMVoxel.inp", "r");
139 if (voxelInpFile == NULL) {
140 printf("neBEMVoxel.inp absent ... assuming OptVoxel = 0 ...\n");
141 OptVoxel = 0;
142 OptStaggerVoxel = 0;
143 } else {
144 fscanf(voxelInpFile, "OptVoxel: %d\n", &OptVoxel);
145 fscanf(voxelInpFile, "OptStaggerVoxel: %d\n", &OptStaggerVoxel);
146 fscanf(voxelInpFile, "Xmin: %le\n", &Voxel.Xmin);
147 fscanf(voxelInpFile, "Xmax: %le\n", &Voxel.Xmax);
148 fscanf(voxelInpFile, "Ymin: %le\n", &Voxel.Ymin);
149 fscanf(voxelInpFile, "Ymax: %le\n", &Voxel.Ymax);
150 fscanf(voxelInpFile, "Zmin: %le\n", &Voxel.Zmin);
151 fscanf(voxelInpFile, "Zmax: %le\n", &Voxel.Zmax);
152 fscanf(voxelInpFile, "XStagger: %le\n", &Voxel.XStagger);
153 fscanf(voxelInpFile, "YStagger: %le\n", &Voxel.YStagger);
154 fscanf(voxelInpFile, "ZStagger: %le\n", &Voxel.ZStagger);
155 fscanf(voxelInpFile, "NbOfXCells: %d\n", &Voxel.NbXCells);
156 fscanf(voxelInpFile, "NbOfYCells: %d\n", &Voxel.NbYCells);
157 fscanf(voxelInpFile, "NbOfZCells: %d\n", &Voxel.NbZCells);
158 fclose(voxelInpFile);
159 } // inputs for Voxel
160
161 // Set up parameters related to 3dMap data export for Garfield++
162 FILE *mapInpFile = fopen("neBEMMap.inp", "r");
163 if (mapInpFile == NULL) {
164 printf("neBEMMap.inp absent ... assuming OptMap = 0 ...\n");
165 OptMap = 0;
166 OptStaggerMap = 0;
167 } else {
168 // While reading the input, OptMap and OptStaggerMap have to be read
169 // first since that will decide whether there is a map and its version.
170 fscanf(mapInpFile, "OptMap: %d\n", &OptMap);
171 fscanf(mapInpFile, "OptStaggerMap: %d\n", &OptStaggerMap);
172 fscanf(mapInpFile, "MapVersion: %9s\n", MapVersion);
173 fscanf(mapInpFile, "Xmin: %le\n", &Map.Xmin);
174 fscanf(mapInpFile, "Xmax: %le\n", &Map.Xmax);
175 fscanf(mapInpFile, "Ymin: %le\n", &Map.Ymin);
176 fscanf(mapInpFile, "Ymax: %le\n", &Map.Ymax);
177 fscanf(mapInpFile, "Zmin: %le\n", &Map.Zmin);
178 fscanf(mapInpFile, "Zmax: %le\n", &Map.Zmax);
179 fscanf(mapInpFile, "XStagger: %le\n", &Map.XStagger);
180 fscanf(mapInpFile, "YStagger: %le\n", &Map.YStagger);
181 fscanf(mapInpFile, "ZStagger: %le\n", &Map.ZStagger);
182 fscanf(mapInpFile, "NbOfXCells: %d\n", &Map.NbXCells);
183 fscanf(mapInpFile, "NbOfYCells: %d\n", &Map.NbYCells);
184 fscanf(mapInpFile, "NbOfZCells: %d\n", &Map.NbZCells);
185 fclose(mapInpFile);
186 } // inputs for 3dMap
187
188 // Set up parameters related to fast volume
189 FILE *fastInpFile = fopen("neBEMFastVol.inp", "r");
190 if (fastInpFile == NULL) {
191 printf("neBEMFastVol.inp absent ... assuming OptFastVol = 0 ...\n");
192 OptFastVol = 0;
194 OptCreateFastPF = 0;
195 OptReadFastPF = 0;
196 FastVol.NbBlocks = 0;
199 } else {
200 fscanf(fastInpFile, "OptFastVol: %d\n", &OptFastVol);
201 fscanf(fastInpFile, "OptStaggerFastVol: %d\n", &OptStaggerFastVol);
202 fscanf(fastInpFile, "OptCreateFastPF: %d\n", &OptCreateFastPF);
203 fscanf(fastInpFile, "OptReadFastPF: %d\n", &OptReadFastPF);
204 fscanf(fastInpFile, "NbPtSkip: %d\n", &NbPtSkip);
205 fscanf(fastInpFile, "NbStgPtSkip: %d\n", &NbStgPtSkip);
206 fscanf(fastInpFile, "LX: %le\n", &FastVol.LX);
207 fscanf(fastInpFile, "LY: %le\n", &FastVol.LY);
208 fscanf(fastInpFile, "LZ: %le\n", &FastVol.LZ);
209 fscanf(fastInpFile, "CornerX: %le\n", &FastVol.CrnrX);
210 fscanf(fastInpFile, "CornerY: %le\n", &FastVol.CrnrY);
211 fscanf(fastInpFile, "CornerZ: %le\n", &FastVol.CrnrZ);
212 fscanf(fastInpFile, "YStagger: %le\n", &FastVol.YStagger);
214 FastVol.YStagger = 0.0; // ignore any non-zero value
215 fscanf(fastInpFile, "NbOfBlocks: %d\n", &FastVol.NbBlocks);
221 for (int block = 1; block <= FastVol.NbBlocks; ++block) {
222 fscanf(fastInpFile, "NbOfXCells: %d\n", &BlkNbXCells[block]);
223 fscanf(fastInpFile, "NbOfYCells: %d\n", &BlkNbYCells[block]);
224 fscanf(fastInpFile, "NbOfZCells: %d\n", &BlkNbZCells[block]);
225 fscanf(fastInpFile, "LZ: %le\n", &BlkLZ[block]);
226 fscanf(fastInpFile, "CornerZ: %le\n", &BlkCrnrZ[block]);
227 } // inputs for blocks
228 fscanf(fastInpFile, "NbOfOmitVols: %d\n", &FastVol.NbOmitVols);
229 if (FastVol.NbOmitVols) {
236 for (int omit = 1; omit <= FastVol.NbOmitVols; ++omit) {
237 fscanf(fastInpFile, "OmitVolLX: %le\n", &OmitVolLX[omit]);
238 fscanf(fastInpFile, "OmitVolLY: %le\n", &OmitVolLY[omit]);
239 fscanf(fastInpFile, "OmitVolLZ: %le\n", &OmitVolLZ[omit]);
240 fscanf(fastInpFile, "OmitVolCornerX: %le\n", &OmitVolCrnrX[omit]);
241 fscanf(fastInpFile, "OmitVolCornerY: %le\n", &OmitVolCrnrY[omit]);
242 fscanf(fastInpFile, "OmitVolCornerZ: %le\n", &OmitVolCrnrZ[omit]);
243 } // inputs for OmitVols
244 } // inputs for OmitVols
245 fscanf(fastInpFile, "NbOfIgnoreVols: %d\n", &FastVol.NbIgnoreVols);
246 if (FastVol.NbIgnoreVols) {
253 for (int ignore = 1; ignore <= FastVol.NbIgnoreVols; ++ignore) {
254 fscanf(fastInpFile, "IgnoreVolLX: %le\n", &IgnoreVolLX[ignore]);
255 fscanf(fastInpFile, "IgnoreVolLY: %le\n", &IgnoreVolLY[ignore]);
256 fscanf(fastInpFile, "IgnoreVolLZ: %le\n", &IgnoreVolLZ[ignore]);
257 fscanf(fastInpFile, "IgnoreVolCornerX: %le\n", &IgnoreVolCrnrX[ignore]);
258 fscanf(fastInpFile, "IgnoreVolCornerY: %le\n", &IgnoreVolCrnrY[ignore]);
259 fscanf(fastInpFile, "IgnoreVolCornerZ: %le\n", &IgnoreVolCrnrZ[ignore]);
260 } // inputs for IgnoreVols
261 } // inputs for IgnoreVols
262 for (int ignore = 1; ignore <= FastVol.NbIgnoreVols; ++ignore) {
263 printf("IgnoreVolLX: %le\n", IgnoreVolLX[ignore]);
264 printf("IgnoreVolLY: %le\n", IgnoreVolLY[ignore]);
265 printf("IgnoreVolLZ: %le\n", IgnoreVolLZ[ignore]);
266 printf("IgnoreVolCornerX: %le\n", IgnoreVolCrnrX[ignore]);
267 printf("IgnoreVolCornerY: %le\n", IgnoreVolCrnrY[ignore]);
268 printf("IgnoreVolCornerZ: %le\n", IgnoreVolCrnrZ[ignore]);
269 } // inputs for IgnoreVols
270 fclose(fastInpFile);
271 } // else fastInpFile
272
273 // Set up parameters related to fixed specification of weighting field
274 FILE *fixedWtInpFile = fopen("neBEMFixedWtField.inp", "r");
275 if (fixedWtInpFile == NULL) {
276 printf(
277 "neBEMFixedWtField.inp absent ... assuming OptFixedWtField = 0 ...\n");
278 OptFixedWtField = 0;
279 FixedWtPotential = 0.0;
280 FixedWtFieldX = 0.0;
281 FixedWtFieldY = 0.0;
282 FixedWtFieldZ = 0.0;
283 } else {
284 fscanf(fixedWtInpFile, "OptFixedWtField: %d\n", &OptFixedWtField);
285 fscanf(fixedWtInpFile, "FixedWtPotential: %lg\n", &FixedWtPotential);
286 fscanf(fixedWtInpFile, "FixedWtFieldX: %lg\n", &FixedWtFieldX);
287 fscanf(fixedWtInpFile, "FixedWtFieldY: %lg\n", &FixedWtFieldY);
288 fscanf(fixedWtInpFile, "FixedWtFieldZ: %lg\n", &FixedWtFieldZ);
289 fclose(fixedWtInpFile);
290 } // else fixedWtInpFile
291
292 // Set up parameters related to weighting field fast volume
293 FILE *fastWtFldInpFile = fopen("neBEMWtFldFastVol.inp", "r");
294 if (fastWtFldInpFile == NULL) {
295 printf(
296 "neBEMWtFldFastVol.inp absent ... assuming OptWtFldFastVol = 0 ...\n");
297 OptWtFldFastVol = 0;
304 } else {
305 fscanf(fastWtFldInpFile, "OptFastVol: %d\n", &OptWtFldFastVol);
306 fscanf(fastWtFldInpFile, "OptStaggerFastVol: %d\n",
308 fscanf(fastWtFldInpFile, "OptCreateFastPF: %d\n", &OptWtFldCreateFastPF);
309 fscanf(fastWtFldInpFile, "OptReadFastPF: %d\n", &OptWtFldReadFastPF);
310 fscanf(fastWtFldInpFile, "NbPtSkip: %d\n", &WtFldNbPtSkip);
311 fscanf(fastWtFldInpFile, "NbStgPtSkip: %d\n", &WtFldNbStgPtSkip);
312 fscanf(fastWtFldInpFile, "LX: %le\n", &WtFldFastVol.LX);
313 fscanf(fastWtFldInpFile, "LY: %le\n", &WtFldFastVol.LY);
314 fscanf(fastWtFldInpFile, "LZ: %le\n", &WtFldFastVol.LZ);
315 fscanf(fastWtFldInpFile, "CornerX: %le\n", &WtFldFastVol.CrnrX);
316 fscanf(fastWtFldInpFile, "CornerY: %le\n", &WtFldFastVol.CrnrY);
317 fscanf(fastWtFldInpFile, "CornerZ: %le\n", &WtFldFastVol.CrnrZ);
318 fscanf(fastWtFldInpFile, "YStagger: %le\n", &WtFldFastVol.YStagger);
320 WtFldFastVol.YStagger = 0.0; // ignore any non-zero value
321 fscanf(fastWtFldInpFile, "NbOfBlocks: %d\n", &WtFldFastVol.NbBlocks);
327 for (int block = 1; block <= WtFldFastVol.NbBlocks; ++block) {
328 fscanf(fastWtFldInpFile, "NbOfXCells: %d\n", &WtFldBlkNbXCells[block]);
329 fscanf(fastWtFldInpFile, "NbOfYCells: %d\n", &WtFldBlkNbYCells[block]);
330 fscanf(fastWtFldInpFile, "NbOfZCells: %d\n", &WtFldBlkNbZCells[block]);
331 fscanf(fastWtFldInpFile, "LZ: %le\n", &WtFldBlkLZ[block]);
332 fscanf(fastWtFldInpFile, "CornerZ: %le\n", &WtFldBlkCrnrZ[block]);
333 } // inputs for blocks
334 fscanf(fastWtFldInpFile, "NbOfOmitVols: %d\n", &WtFldFastVol.NbOmitVols);
342 for (int omit = 1; omit <= WtFldFastVol.NbOmitVols; ++omit) {
343 fscanf(fastWtFldInpFile, "OmitVolLX: %le\n", &WtFldOmitVolLX[omit]);
344 fscanf(fastWtFldInpFile, "OmitVolLY: %le\n", &WtFldOmitVolLY[omit]);
345 fscanf(fastWtFldInpFile, "OmitVolLZ: %le\n", &WtFldOmitVolLZ[omit]);
346 fscanf(fastWtFldInpFile, "OmitVolCornerX: %le\n",
347 &WtFldOmitVolCrnrX[omit]);
348 fscanf(fastWtFldInpFile, "OmitVolCornerY: %le\n",
349 &WtFldOmitVolCrnrY[omit]);
350 fscanf(fastWtFldInpFile, "OmitVolCornerZ: %le\n",
351 &WtFldOmitVolCrnrZ[omit]);
352 } // inputs for OmitVols
353 } // inputs for OmitVols
354 fscanf(fastWtFldInpFile, "NbOfIgnoreVols: %d\n",
363 for (int ignore = 1; ignore <= WtFldFastVol.NbIgnoreVols; ++ignore) {
364 fscanf(fastWtFldInpFile, "IgnoreVolLX: %le\n",
365 &WtFldIgnoreVolLX[ignore]);
366 fscanf(fastWtFldInpFile, "IgnoreVolLY: %le\n",
367 &WtFldIgnoreVolLY[ignore]);
368 fscanf(fastWtFldInpFile, "IgnoreVolLZ: %le\n",
369 &WtFldIgnoreVolLZ[ignore]);
370 fscanf(fastWtFldInpFile, "IgnoreVolCornerX: %le\n",
371 &WtFldIgnoreVolCrnrX[ignore]);
372 fscanf(fastWtFldInpFile, "IgnoreVolCornerY: %le\n",
373 &WtFldIgnoreVolCrnrY[ignore]);
374 fscanf(fastWtFldInpFile, "IgnoreVolCornerZ: %le\n",
375 &WtFldIgnoreVolCrnrZ[ignore]);
376 } // inputs for IgnoreVols
377 } // inputs for IgnoreVols
378 for (int ignore = 1; ignore <= WtFldFastVol.NbIgnoreVols; ++ignore) {
379 printf("WtFldIgnoreVolLX: %le\n", WtFldIgnoreVolLX[ignore]);
380 printf("WtFldIgnoreVolLY: %le\n", WtFldIgnoreVolLY[ignore]);
381 printf("WtFldIgnoreVolLZ: %le\n", WtFldIgnoreVolLZ[ignore]);
382 printf("WtFldIgnoreVolCornerX: %le\n", WtFldIgnoreVolCrnrX[ignore]);
383 printf("WtFldIgnoreVolCornerY: %le\n", WtFldIgnoreVolCrnrY[ignore]);
384 printf("WtFldIgnoreVolCornerZ: %le\n", WtFldIgnoreVolCrnrZ[ignore]);
385 } // inputs for IgnoreVols
386 fclose(fastWtFldInpFile);
387 } // else fastWtFldInpFile
388
389 printf("neBEM initialized ...\n");
390 fflush(stdout);
391 sleep(3); // wait for three seconds so that the user gets time to react
392
393 neBEMState = 1; // state 1 implied initialization of neBEM completed
394
395 // announce success - later, add the name of the calling code
396 return (0);
397} // neBEMInitialize ends
ISLESGLOBAL char ISLESVersion[10]
Definition: Isles.h:30
int CreateDirStr(void)
INTFACEGLOBAL int NbThreads
INTFACEGLOBAL char DeviceInputFile[256]
INTFACEGLOBAL int neBEMGetInputsFromFiles(void)
INTFACEGLOBAL int OptDeviceFile
INTFACEGLOBAL int neBEMSetDefaults(void)
neBEMGLOBAL int OptWtFldReadFastPF
Definition: neBEM.h:475
neBEMGLOBAL int OptStaggerFastVol
Definition: neBEM.h:404
neBEMGLOBAL int OptVoxel
Definition: neBEM.h:354
neBEMGLOBAL double * OmitVolCrnrZ
Definition: neBEM.h:434
neBEMGLOBAL int * WtFldBlkNbYCells
Definition: neBEM.h:494
neBEMGLOBAL double * WtFldOmitVolLX
Definition: neBEM.h:498
neBEMGLOBAL double FixedWtFieldZ
Definition: neBEM.h:469
neBEMGLOBAL int OptWtFldFastVol
Definition: neBEM.h:472
neBEMGLOBAL MapVol Map
Definition: neBEM.h:396
neBEMGLOBAL double * OmitVolLX
Definition: neBEM.h:429
neBEMGLOBAL int OptStaggerMap
Definition: neBEM.h:379
neBEMGLOBAL char neBEMVersion[10]
Definition: neBEM.h:34
neBEMGLOBAL double * WtFldOmitVolCrnrZ
Definition: neBEM.h:503
neBEMGLOBAL int * WtFldBlkNbXCells
Definition: neBEM.h:493
neBEMGLOBAL double * WtFldIgnoreVolCrnrX
Definition: neBEM.h:507
neBEMGLOBAL double * BlkCrnrZ
Definition: neBEM.h:428
neBEMGLOBAL int OptCreateFastPF
Definition: neBEM.h:405
neBEMGLOBAL double * WtFldOmitVolCrnrY
Definition: neBEM.h:502
neBEMGLOBAL int OptWtFldStaggerFastVol
Definition: neBEM.h:473
neBEMGLOBAL FastAlgoVol FastVol
Definition: neBEM.h:422
neBEMGLOBAL double * WtFldOmitVolLZ
Definition: neBEM.h:500
neBEMGLOBAL int WtFldNbPtSkip
Definition: neBEM.h:476
neBEMGLOBAL int PrimAfter
Definition: neBEM.h:351
neBEMGLOBAL double * BlkLZ
Definition: neBEM.h:427
neBEMGLOBAL double * OmitVolCrnrX
Definition: neBEM.h:432
neBEMGLOBAL double * OmitVolLZ
Definition: neBEM.h:431
neBEMGLOBAL int NbStgPtSkip
Definition: neBEM.h:408
neBEMGLOBAL double * IgnoreVolCrnrY
Definition: neBEM.h:439
neBEMGLOBAL double FixedWtFieldX
Definition: neBEM.h:467
neBEMGLOBAL int OptWtFldCreateFastPF
Definition: neBEM.h:474
neBEMGLOBAL char MapVersion[10]
Definition: neBEM.h:380
neBEMGLOBAL int NbPtSkip
Definition: neBEM.h:407
neBEMGLOBAL double * WtFldIgnoreVolLZ
Definition: neBEM.h:506
neBEMGLOBAL int OptMap
Definition: neBEM.h:378
neBEMGLOBAL WtFldFastAlgoVol WtFldFastVol
Definition: neBEM.h:491
neBEMGLOBAL double * OmitVolCrnrY
Definition: neBEM.h:433
neBEMGLOBAL double * WtFldBlkLZ
Definition: neBEM.h:496
neBEMGLOBAL double * WtFldBlkCrnrZ
Definition: neBEM.h:497
neBEMGLOBAL double * OmitVolLY
Definition: neBEM.h:430
neBEMGLOBAL int * BlkNbXCells
Definition: neBEM.h:424
neBEMGLOBAL int OptReadFastPF
Definition: neBEM.h:406
neBEMGLOBAL int * BlkNbZCells
Definition: neBEM.h:426
neBEMGLOBAL double * WtFldOmitVolLY
Definition: neBEM.h:499
neBEMGLOBAL int * WtFldBlkNbZCells
Definition: neBEM.h:495
neBEMGLOBAL double FixedWtFieldY
Definition: neBEM.h:468
neBEMGLOBAL int WtFldNbStgPtSkip
Definition: neBEM.h:477
neBEMGLOBAL double * IgnoreVolLY
Definition: neBEM.h:436
neBEMGLOBAL double FixedWtPotential
Definition: neBEM.h:466
neBEMGLOBAL double * WtFldIgnoreVolCrnrY
Definition: neBEM.h:508
neBEMGLOBAL double * IgnoreVolLX
Definition: neBEM.h:435
neBEMGLOBAL double * WtFldIgnoreVolCrnrZ
Definition: neBEM.h:509
neBEMGLOBAL double * WtFldIgnoreVolLY
Definition: neBEM.h:505
neBEMGLOBAL double LengthScale
Definition: neBEM.h:238
neBEMGLOBAL int OptFixedWtField
Definition: neBEM.h:465
neBEMGLOBAL double * IgnoreVolLZ
Definition: neBEM.h:437
neBEMGLOBAL double * WtFldOmitVolCrnrX
Definition: neBEM.h:501
neBEMGLOBAL VoxelVol Voxel
Definition: neBEM.h:371
neBEMGLOBAL double * WtFldIgnoreVolLX
Definition: neBEM.h:504
neBEMGLOBAL double * IgnoreVolCrnrX
Definition: neBEM.h:438
neBEMGLOBAL double * IgnoreVolCrnrZ
Definition: neBEM.h:440
neBEMGLOBAL int * BlkNbYCells
Definition: neBEM.h:425
neBEMGLOBAL int OptStaggerVoxel
Definition: neBEM.h:355
int * ivector(long nl, long nh)
Definition: nrutil.c:32
double * dvector(long nl, long nh)
Definition: nrutil.c:63
double CrnrX
Definition: neBEM.h:414
int NbBlocks
Definition: neBEM.h:418
double CrnrZ
Definition: neBEM.h:416
double CrnrY
Definition: neBEM.h:415
double LZ
Definition: neBEM.h:413
int NbOmitVols
Definition: neBEM.h:419
int NbIgnoreVols
Definition: neBEM.h:420
double YStagger
Definition: neBEM.h:417
double LY
Definition: neBEM.h:412
double LX
Definition: neBEM.h:411
int NbYCells
Definition: neBEM.h:393
double Ymax
Definition: neBEM.h:386
double YStagger
Definition: neBEM.h:390
double Zmax
Definition: neBEM.h:388
double ZStagger
Definition: neBEM.h:391
int NbZCells
Definition: neBEM.h:394
double Xmax
Definition: neBEM.h:384
double Xmin
Definition: neBEM.h:383
double XStagger
Definition: neBEM.h:389
double Zmin
Definition: neBEM.h:387
int NbXCells
Definition: neBEM.h:392
double Ymin
Definition: neBEM.h:385
double Zmax
Definition: neBEM.h:363
double YStagger
Definition: neBEM.h:365
double Xmax
Definition: neBEM.h:359
int NbXCells
Definition: neBEM.h:367
double Ymax
Definition: neBEM.h:361
double Zmin
Definition: neBEM.h:362
double ZStagger
Definition: neBEM.h:366
double XStagger
Definition: neBEM.h:364
int NbZCells
Definition: neBEM.h:369
int NbYCells
Definition: neBEM.h:368
double Xmin
Definition: neBEM.h:358
double Ymin
Definition: neBEM.h:360
int NbIgnoreVols
Definition: neBEM.h:489
int NbOmitVols
Definition: neBEM.h:488
double LZ
Definition: neBEM.h:482
double YStagger
Definition: neBEM.h:486
double LY
Definition: neBEM.h:481
double LX
Definition: neBEM.h:480
double CrnrY
Definition: neBEM.h:484
double CrnrX
Definition: neBEM.h:483
double CrnrZ
Definition: neBEM.h:485

◆ neBEMKnownCharges()

INTFACEGLOBAL int neBEMKnownCharges ( void  )

Definition at line 2244 of file neBEMInterface.c.

2244 {
2245 int debugFn = 0;
2246
2247 /*
2248 // How to check that the pointer points to a valid function?
2249 if(Pt2UserFunction == NULL)
2250 {
2251 printf("Not a valid function ... returning ...\n");
2252 return(-1);
2253 }
2254 else
2255 {
2256 // printf("Pt2UserFunction points to %p\n", Pt2UserFunction);
2257 }
2258
2259 // Status of the known charges conditions is meaningful only after the
2260 // element discretization has been completed, i.e., beyond the 5th state.
2261 if(neBEMState >= 5)
2262 { // the following function is declared in the Interface.h.
2263 int fstatus = (*Pt2UserFunction)(); // user to supply function
2264 if(fstatus != 0)
2265 {
2266 neBEMMessage("neBEMKnownCharges - Pt2UserFunction");
2267 return -1;
2268 }
2269 if(neBEMState > 5) // assume LHS and inversion to be over?
2270 neBEMState = 8;
2271 }
2272 else
2273 {
2274 printf("Known charges are meaningful only beyond state 4 ...\n");
2275 printf("returning ...\n");
2276 return(-1);
2277 }
2278 */
2279
2280 // Set up parameters related to known charge calculations
2281 // Electrons and ions can be distributed as points, lines, areas or volumes.
2282 {
2283 FILE *KnChInpFile = fopen("neBEMKnCh.inp", "r");
2284 if (KnChInpFile == NULL) {
2285 OptKnCh = 0;
2286 printf(
2287 "neBEMKnCh.inp absent ... assuming absence of known charges ...\n");
2288 } else {
2289 fscanf(KnChInpFile, "OptKnCh: %d\n", &OptKnCh);
2290 if (1) printf("OptKnCh: %d\n", OptKnCh);
2291
2292 if (!OptKnCh) printf("OptKnCh = 0 ... assuming no known charges ...\n");
2293
2294 if (OptKnCh) {
2295 char PointKnChFile[256];
2296 char LineKnChFile[256];
2297 char AreaKnChFile[256];
2298 char VolumeKnChFile[256];
2299 double KnChFactor;
2300
2301 fscanf(KnChInpFile, "NbPointsKnCh: %d\n", &NbPointsKnCh);
2302 fscanf(KnChInpFile, "PointKnChFile: %255s\n", PointKnChFile);
2303 fscanf(KnChInpFile, "NbLinesKnCh: %d\n", &NbLinesKnCh);
2304 fscanf(KnChInpFile, "LineKnChFile: %255s\n", LineKnChFile);
2305 fscanf(KnChInpFile, "NbAreasKnCh: %d\n", &NbAreasKnCh);
2306 fscanf(KnChInpFile, "AreaKnChFile: %255s\n", AreaKnChFile);
2307 fscanf(KnChInpFile, "NbVolumesKnCh: %d\n", &NbVolumesKnCh);
2308 fscanf(KnChInpFile, "VolumeKnChFile: %255s\n", VolumeKnChFile);
2309 fscanf(KnChInpFile, "KnChFactor: %lg\n", &KnChFactor);
2310 if (1) {
2311 printf("NbPointsKnCh: %d\n", NbPointsKnCh);
2312 printf("PointKnChFile: %s\n", PointKnChFile);
2313 printf("NbLinesKnCh: %d\n", NbLinesKnCh);
2314 printf("LineKnChFile: %s\n", LineKnChFile);
2315 printf("NbAreasKnCh: %d\n", NbAreasKnCh);
2316 printf("AreaKnChFile: %s\n", AreaKnChFile);
2317 printf("NbVolumesKnCh: %d\n", NbVolumesKnCh);
2318 printf("VolumeKnChFile: %s\n", VolumeKnChFile);
2319 printf("KnChFactor: %lg\n", KnChFactor);
2320 }
2321
2322 if (NbPointsKnCh) { // Points
2323 if (debugFn)
2324 printf("No. of points with known charges: %d\n", NbPointsKnCh);
2325
2326 FILE *fptrPointKnChFile = fopen(PointKnChFile, "r");
2327 if (fptrPointKnChFile == NULL) {
2328 neBEMMessage("PointKnCh file absent ... returning\n");
2329 return -10;
2330 }
2331
2332 PointKnChArr =
2333 (PointKnCh *)malloc((NbPointsKnCh + 1) * sizeof(PointKnCh));
2334 if (PointKnChArr == NULL) {
2335 neBEMMessage("Memory allocation failed ... returning\n");
2336 fclose(fptrPointKnChFile);
2337 return -10;
2338 }
2339 for (int point = 0; point <= NbPointsKnCh; ++point) {
2340 // CHECK!!! ele limits start from 0, but all else from 1 to ...
2341 PointKnChArr[point].Nb = 0;
2342 PointKnChArr[point].P.X = 0.0;
2343 PointKnChArr[point].P.Y = 0.0;
2344 PointKnChArr[point].P.Z = 0.0;
2345 PointKnChArr[point].Assigned = 0.0;
2346 }
2347
2348 char header[256];
2349 fgets(header, 256, fptrPointKnChFile); // header
2350
2351 for (int point = 1; point <= NbPointsKnCh; ++point) {
2352 fscanf(fptrPointKnChFile, "%d %lg %lg %lg %lg\n",
2353 &PointKnChArr[point].Nb, &PointKnChArr[point].P.X,
2354 &PointKnChArr[point].P.Y, &PointKnChArr[point].P.Z,
2355 &PointKnChArr[point].Assigned);
2356
2357 PointKnChArr[point].Assigned *= KnChFactor;
2358 }
2359
2360 fclose(fptrPointKnChFile);
2361 if (debugFn) printf("done for all points\n");
2362 } // if NbPointsKnCh: Inputs and calculations for points ends
2363
2364 if (NbLinesKnCh) { // Lines
2365 if (debugFn)
2366 printf("No. of lines with known charges: %d\n", NbLinesKnCh);
2367
2368 FILE *fptrLineKnChFile = fopen(LineKnChFile, "r");
2369 if (fptrLineKnChFile == NULL) {
2370 neBEMMessage("LineKnCh file absent ... returning\n");
2371 return -10;
2372 }
2373
2374 LineKnChArr =
2375 (LineKnCh *)malloc((NbLinesKnCh + 1) * sizeof(LineKnCh));
2376 if (LineKnChArr == NULL) {
2377 neBEMMessage("Memory allocation failed ... returning\n");
2378 fclose(fptrLineKnChFile);
2379 return -10;
2380 }
2381 for (int line = 0; line <= NbLinesKnCh; ++line) {
2382 // CHECK!!! ele limits start from 0, but all else from 1 to ...
2383 LineKnChArr[line].Nb = 0;
2384 LineKnChArr[line].Start.X = 0.0;
2385 LineKnChArr[line].Start.Y = 0.0;
2386 LineKnChArr[line].Start.Z = 0.0;
2387 LineKnChArr[line].Stop.X = 0.0;
2388 LineKnChArr[line].Stop.Y = 0.0;
2389 LineKnChArr[line].Stop.Z = 0.0;
2390 LineKnChArr[line].Radius = 0.0;
2391 LineKnChArr[line].Assigned = 0.0;
2392 }
2393
2394 char header[256];
2395 fgets(header, 256, fptrLineKnChFile); // header
2396
2397 for (int line = 1; line <= NbLinesKnCh; ++line) {
2398 fscanf(fptrLineKnChFile, "%d %lg %lg %lg %lg %lg %lg %lg %lg\n",
2399 &LineKnChArr[line].Nb, &LineKnChArr[line].Start.X,
2400 &LineKnChArr[line].Start.Y, &LineKnChArr[line].Start.Z,
2401 &LineKnChArr[line].Stop.X, &LineKnChArr[line].Stop.Y,
2402 &LineKnChArr[line].Stop.Z, &LineKnChArr[line].Radius,
2403 &LineKnChArr[line].Assigned);
2404 LineKnChArr[line].Assigned *= KnChFactor;
2405 }
2406
2407 fclose(fptrLineKnChFile);
2408 if (debugFn) printf("done for all lines\n");
2409 } // if NbLinesKnCh: Inputs and calculations for lines ends
2410
2411 if (NbAreasKnCh) { // Areas
2412 if (debugFn)
2413 printf("No. of areas with known charges: %d\n", NbAreasKnCh);
2414
2415 FILE *fptrAreaKnChFile = fopen(AreaKnChFile, "r");
2416 if (fptrAreaKnChFile == NULL) {
2417 neBEMMessage("AreaKnCh file absent ... returning\n");
2418 return -10;
2419 }
2420
2421 AreaKnChArr =
2422 (AreaKnCh *)malloc((NbAreasKnCh + 1) * sizeof(AreaKnCh));
2423 if (AreaKnChArr == NULL) {
2424 neBEMMessage("Memory allocation failed ... returning\n");
2425 fclose(fptrAreaKnChFile);
2426 return -10;
2427 }
2428 for (int area = 0; area <= NbAreasKnCh; ++area) {
2429 // CHECK!!! ele limits start from 0, but all else from 1 to ...
2430 AreaKnChArr[area].Nb = 0;
2431 AreaKnChArr[area].NbVertices = 0;
2432 for (int vert = 1; vert <= (AreaKnChArr + area - 1)->NbVertices;
2433 ++vert) {
2434 (AreaKnChArr + area - 1)->Vertex[vert].X = 0.0;
2435 (AreaKnChArr + area - 1)->Vertex[vert].Y = 0.0;
2436 (AreaKnChArr + area - 1)->Vertex[vert].Z = 0.0;
2437 }
2438 AreaKnChArr[area].Assigned = 0.0;
2439 }
2440
2441 char header[256];
2442 fgets(header, 256, fptrAreaKnChFile); // header
2443
2444 for (int area = 1; area <= NbAreasKnCh; ++area) {
2445 fscanf(fptrAreaKnChFile, "%d %d %le\n",
2446 &(AreaKnChArr + area - 1)->Nb,
2447 &(AreaKnChArr + area - 1)->NbVertices,
2448 &(AreaKnChArr + area - 1)->Assigned);
2449 for (int vert = 1; vert <= (AreaKnChArr + area - 1)->NbVertices;
2450 ++vert) {
2451 fscanf(fptrAreaKnChFile, "%le %le %le\n",
2452 &(AreaKnChArr + area - 1)->Vertex[vert].X,
2453 &(AreaKnChArr + area - 1)->Vertex[vert].Y,
2454 &(AreaKnChArr + area - 1)->Vertex[vert].Z);
2455 }
2456 AreaKnChArr[area].Assigned *= KnChFactor;
2457 }
2458
2459 fclose(fptrAreaKnChFile);
2460 if (debugFn) printf("done for all areas\n");
2461 } // if AreasKnCh: Inputs and calculations for areas ends
2462
2463 if (NbVolumesKnCh) { // Volumes
2464 if (debugFn)
2465 printf("No. of volumes with known charges: %d\n", NbVolumesKnCh);
2466
2467 FILE *fptrVolumeKnChFile = fopen(VolumeKnChFile, "r");
2468 if (fptrVolumeKnChFile == NULL) {
2469 neBEMMessage("VolumeKnCh file absent ... returning\n");
2470 return -10;
2471 }
2472
2474 (VolumeKnCh *)malloc((NbVolumesKnCh + 1) * sizeof(VolumeKnCh));
2475 if (VolumeKnChArr == NULL) {
2476 neBEMMessage("Memory allocation failed ... returning\n");
2477 fclose(fptrVolumeKnChFile);
2478 return -10;
2479 }
2480 for (int volume = 0; volume <= NbVolumesKnCh; ++volume) {
2481 // CHECK!!! ele limits start from 0, but all else from 1 to ...
2482 VolumeKnChArr[volume].Nb = 0;
2483 VolumeKnChArr[volume].NbVertices = 0;
2484 for (int vert = 1; vert <= (VolumeKnChArr + volume - 1)->NbVertices;
2485 ++vert) {
2486 (VolumeKnChArr + volume - 1)->Vertex[vert].X = 0.0;
2487 (VolumeKnChArr + volume - 1)->Vertex[vert].Y = 0.0;
2488 (VolumeKnChArr + volume - 1)->Vertex[vert].Z = 0.0;
2489 }
2490 VolumeKnChArr[volume].Assigned = 0.0;
2491 }
2492
2493 char header[256];
2494 fgets(header, 256, fptrVolumeKnChFile); // header
2495
2496 for (int volume = 1; volume <= NbVolumesKnCh; ++volume) {
2497 fscanf(fptrVolumeKnChFile, "%d %d %le\n",
2498 &(VolumeKnChArr + volume - 1)->Nb,
2499 &(VolumeKnChArr + volume - 1)->NbVertices,
2500 &(VolumeKnChArr + volume - 1)->Assigned);
2501 for (int vert = 1; vert <= (VolumeKnChArr + volume - 1)->NbVertices;
2502 ++vert) {
2503 fscanf(fptrVolumeKnChFile, "%le %le %le\n",
2504 &(VolumeKnChArr + volume - 1)->Vertex[vert].X,
2505 &(VolumeKnChArr + volume - 1)->Vertex[vert].Y,
2506 &(VolumeKnChArr + volume - 1)->Vertex[vert].Z);
2507 }
2508 VolumeKnChArr[volume].Assigned *= KnChFactor;
2509 }
2510
2511 fclose(fptrVolumeKnChFile);
2512 if (debugFn) printf("done for all volumes\n");
2513 } // if NbVolumesKnCh: Inputs and calculations for volumes ends
2514
2515 } // if OptKnCh
2516
2517 fclose(KnChInpFile);
2518 } // else KnChInpFile
2519 } // parameters related to known charge calculations
2520
2521 return (0);
2522} // neBEMKnownCharges ends
neBEMGLOBAL PointKnCh * PointKnChArr
Definition: neBEM.h:181
neBEMGLOBAL VolumeKnCh * VolumeKnChArr
Definition: neBEM.h:216
neBEMGLOBAL int NbPointsKnCh
Definition: neBEM.h:172
neBEMGLOBAL AreaKnCh * AreaKnChArr
Definition: neBEM.h:204
neBEMGLOBAL int NbLinesKnCh
Definition: neBEM.h:172
neBEMGLOBAL int NbVolumesKnCh
Definition: neBEM.h:172
neBEMGLOBAL int NbAreasKnCh
Definition: neBEM.h:172
neBEMGLOBAL LineKnCh * LineKnChArr
Definition: neBEM.h:192
neBEMGLOBAL int OptKnCh
Definition: neBEM.h:52
double Assigned
Definition: neBEM.h:201
int NbVertices
Definition: neBEM.h:199
int Nb
Definition: neBEM.h:198
Point3D Stop
Definition: neBEM.h:187
int Nb
Definition: neBEM.h:185
double Radius
Definition: neBEM.h:188
Point3D Start
Definition: neBEM.h:186
double Assigned
Definition: neBEM.h:189
Point3D P
Definition: neBEM.h:177
double Assigned
Definition: neBEM.h:178
int Nb
Definition: neBEM.h:176
int Nb
Definition: neBEM.h:210
int NbVertices
Definition: neBEM.h:211
double Assigned
Definition: neBEM.h:213

Referenced by ComputeSolution().

◆ neBEMMessage()

◆ neBEMPrepareWeightingField()

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

Definition at line 4605 of file neBEMInterface.c.

4605 {
4606 static int IdWtField = 0;
4607
4608 if (neBEMState < 7) {
4609 printf(
4610 "neBEMPrepareWeightingField: Weighting computations only meaningful "
4611 "beyond neBEMState 7 ...\n");
4612 return -1;
4613 }
4614
4615 // Find first free slot
4616 const int MaxWtField = 100; // used again while deallocating these memories
4617 if (WtFieldChDen == NULL)
4618 WtFieldChDen = (double **)malloc(MaxWtField * sizeof(double *));
4619 if (AvWtChDen == NULL)
4620 AvWtChDen = (double **)malloc(MaxWtField * sizeof(double *));
4621
4622 ++IdWtField;
4623 if (IdWtField >= MaxWtField) {
4624 printf(
4625 "neBEMPrepareWeightingField: reached MaxWtField weighting fields.\n");
4626 return -1;
4627 }
4628
4629 // Allocate a new column to store this solution set
4630 WtFieldChDen[IdWtField] = (double *)malloc((NbElements + 2) * sizeof(double));
4631 AvWtChDen[IdWtField] = (double *)malloc((NbPrimitives + 2) * sizeof(double));
4632
4633 int fstatus =
4634 WeightingFieldSolution(nprim, primlist, WtFieldChDen[IdWtField]);
4635 if (!fstatus) // estimate primitive related average wt field charge densities
4636 {
4637 // OMPCheck - may be parallelized
4638 for (int prim = 1; prim <= NbPrimitives; ++prim) {
4639 double area = 0.0; // need area of the primitive as well!
4640 AvWtChDen[IdWtField][prim] = 0.0;
4641
4642 for (int ele = ElementBgn[prim]; ele <= ElementEnd[prim]; ++ele) {
4643 area += (EleArr + ele - 1)->G.dA;
4644 AvWtChDen[IdWtField][prim] +=
4645 WtFieldChDen[IdWtField][ele] * (EleArr + ele - 1)->G.dA;
4646 }
4647
4648 AvWtChDen[IdWtField][prim] /= area;
4649 }
4650 }
4651 if (fstatus != 0) {
4652 neBEMMessage("neBEMPrepareWeightingField - WeightingFieldSolution");
4653 return -1;
4654 }
4655
4656 return IdWtField;
4657} // neBEMPrepareWeightingField ends
int WeightingFieldSolution(int NbPrimsWtField, int PrimListWtField[], double solnarray[])
Definition: neBEM.c:3688

◆ neBEMReadGeometry()

INTFACEGLOBAL int neBEMReadGeometry ( void  )

Definition at line 406 of file neBEMInterface.c.

406 {
407 int dbgFn = 0;
408 int fstatus;
409
410 startClock = clock();
411
412 // For a model that was defined before and for which data was stored in a file
413 if ((!NewModel) && (!NewBC) && (OptStorePrimitives)) {
414 fstatus = ReadPrimitives();
415 if (fstatus) {
416 neBEMMessage("neBEMReadGeometry - problem reading stored Primitives.\n");
417 return -1;
418 }
419 neBEMState = 3; // primitives read in after initialization and Nbs
420 return 0;
421 }
422
423 printf("geometry inputs ...\n");
424 if (neBEMState != 1) {
425 printf("reading geometry possible only after initialization ...\n");
426 return -1;
427 }
428
431 if (NbPrimitives == 0) {
432 // nothing to do - return control to calling routine
433 neBEMMessage("neBEMReadGeometry - no primitive.\n");
434 return (-1); // for the time being
435 }
436
437 NbSurfs = 0;
438 NbWires = 0;
439 neBEMState = 2;
440
441 // Allocate memory for storing the geometry primitives till the elements are
442 // created.
443 // Explicit storage of these variables related to primitives may not be
444 // necessary if the elements are created immediately after a primitive is read
445 // in.
446 // MaxNbVertices = 4; // value specified through SetDefaults or init files
447
448 // neBEM has been initialized, NbPrimitives set
451 OrgnlToEffPrim = imatrix(1, NbPrimitives, 0, 2); // 0 init, 1 intrfc, 2 rmv
460 Radius = dvector(1, NbPrimitives); // can lead to a little memory misuse
464 PrimDC = (DirnCosn3D *)malloc(NbPrimitives * sizeof(DirnCosn3D));
469 NbWireSeg = ivector(1, NbPrimitives); // little memory misuse
514
515 // Loop over the primitives - major loop
516 int nvertex, volref1, volref2, volmax = 0;
517#ifdef __cplusplus
518 std::vector<double> xvert(MaxNbVertices, 0.);
519 std::vector<double> yvert(MaxNbVertices, 0.);
520 std::vector<double> zvert(MaxNbVertices, 0.);
521#else
522 double xvert[MaxNbVertices], yvert[MaxNbVertices], zvert[MaxNbVertices];
523#endif
524 double xnorm, ynorm, znorm; // in case of wire , radius is read as xnorm
525 for (int prim = 1; prim <= NbPrimitives; ++prim) {
526#ifdef __cplusplus
527 fstatus = neBEMGetPrimitive(prim, &nvertex,
528 xvert.data(), yvert.data(), zvert.data(),
529 &xnorm, &ynorm, &znorm, &volref1, &volref2);
530#else
531 fstatus = neBEMGetPrimitive(prim, &nvertex, xvert, yvert, zvert, // arrays
532 &xnorm, &ynorm, &znorm, &volref1, &volref2);
533#endif
534 if (fstatus != 0) {
535 neBEMMessage("neBEMReadGeometry - neBEMGetPrimitve");
536 return -1;
537 }
538 if (volmax < volref1) {
539 volmax = volref1;
540 } // maxm nb of volumes
541 if (volmax < volref2) {
542 volmax = volref2;
543 } // maxm nb of volumes
544
545 if (nvertex > MaxNbVertices) {
546 printf("Number of vertices for primitive %d exceeds %d!\n", prim,
548 printf("Returning to garfield ...\n");
549 return (-1);
550 }
551
552 PrimType[prim] = nvertex; // wire:2, triangle:3, rectangle:4
553 NbVertices[prim] = nvertex;
554 for (int vert = 0; vert < NbVertices[prim]; ++vert) {
555 XVertex[prim][vert] = xvert[vert];
556 YVertex[prim][vert] = yvert[vert];
557 ZVertex[prim][vert] = zvert[vert];
558 }
559 if (PrimType[prim] == 2) {
560 // wire
561 XNorm[prim] = 0.0; // modulus not 1 - an absurd trio!
562 YNorm[prim] = 0.0;
563 ZNorm[prim] = 0.0;
564 Radius[prim] = xnorm;
565 }
566 if ((PrimType[prim] == 3) || (PrimType[prim] == 4)) {
567 XNorm[prim] = xnorm;
568 YNorm[prim] = ynorm;
569 ZNorm[prim] = znorm;
570 Radius[prim] = 0.0; // absurd radius!
571 }
572 VolRef1[prim] = volref1;
573 VolRef2[prim] = volref2;
574
575 // feedback for user begins - suppress later
577 printf("neBEM:\tprimitive %d between volumes %d, %d has %d vertices\n",
578 prim, volref1, volref2, nvertex);
579 for (int ivertex = 0; ivertex < nvertex; ivertex++) {
580 printf("\tnode %d (%g,%g,%g)\n", ivertex, xvert[ivertex],
581 yvert[ivertex], zvert[ivertex]);
582 }
583 printf("\tnormal vector: (%g,%g,%g)\n", xnorm, ynorm, znorm);
584 }
585 // feedback for user ends - suppress later
586
587 // Now look for the volume related information for this primitve
588 // This is obtained from the specified volume references
589 // volref1 refers to the volume itself
590 // volref2 describes the volume in the direction of the +ve normal
591 // Note that materials from 1 to 10 are conductors and
592 // from 11 to 20 are
593 // dielectrics
594 int shape1, material1, boundarytype1;
595 double epsilon1, potential1, charge1;
596 if (volref1 == -1) {
597 // Must be an error, since no device is made of vacuum
598 neBEMMessage("neBEMReadGeometry - volref1 = -1!");
599 return -1;
600 } else {
601 neBEMVolumeDescription(volref1, &shape1, &material1, &epsilon1,
602 &potential1, &charge1, &boundarytype1);
603 }
605 printf("\tvolref1: %d\n", volref1);
606 printf("\t\tboundarytype1: %d, shape1: %d, material1: %d\n",
607 boundarytype1, shape1, material1);
608 printf("\t\tepsilon1: %lg, potential1: %lg, charge1: %lg\n", epsilon1,
609 potential1, charge1);
610 }
611 // in the -ve normal direction - properties of the external volume
612 int shape2, material2, boundarytype2;
613 double epsilon2, potential2, charge2;
614 if (volref2 == -1) {
615 shape2 = 0;
616 material2 = 11;
617 epsilon2 = 1.0;
618 potential2 = 0.0;
619 charge2 = 0.0;
620 boundarytype2 = 0;
621 } else {
622 neBEMVolumeDescription(volref2, &shape2, &material2, &epsilon2,
623 &potential2, &charge2, &boundarytype2);
624 }
626 printf("\tvolref2: %d\n", volref2);
627 printf("\t\tboundarytype2: %d, shape2: %d, material2: %d\n",
628 boundarytype2, shape2, material2);
629 printf("\t\tepsilon2: %lg, potential2: %lg, charge2: %lg\n", epsilon2,
630 potential2, charge2);
631 }
632
633 // Put default values to variables that depend on the interface type
634 // Is there any risk involved in putting these defaults?
635 // At present, they even seem necessary. For example, for floating
636 // conductors or dielectric-dielectric interface, the formulation requires
637 // that the RHS is zero (may be modified by the effect of known charges).
638 Epsilon1[prim] = epsilon1;
639 Epsilon2[prim] = epsilon2; // 1: self, 2: external
640 ApplPot[prim] = 0.0;
641 Lambda[prim] = 0.0;
642 ApplCh[prim] = 0.0;
643
644 // BoundaryTypes:
645 // --------------
646 // Vacuum: 0
647 // Conductor at specified potential: 1
648 // Conductor with a specified charge: 2
649 // Floating conductor (zero charge, perpendicular E): 3
650 // Dielectric intrface (plastic-plastic) without "manual" charge: 4
651 // Dielectric with surface charge (plastic-gas, typically): 5
652 // Symmetry boundary, E parallel: 6 (may not be necessary)
653 // Symmetry boundary, E perpendicular: 7 (may not be necessary)
654
655 // InterfaceType:
656 // --------------
657 // To be skipped: 0
658 // Conductor-dielectric: 1
659 // Conductor with known charge: 2
660 // Conductor at floating potential: 3
661 // Dielectric-dielectric: 4
662 // Dielectric with given surface charge: 5
663 // Check dielectric-dielectric formulation in
664 // (NumSolnOfBIEforMolES_JPBardhan.pdf):
665 // Numerical solution of boundary-integral equations for molecular
666 // electrostatics,
667 // by Jaydeep P. Bardhan,
668 // THE JOURNAL OF CHEMICAL PHYSICS 130, 094102 (2009)
669
670 switch (boundarytype1) { // the volume itself is volref1
671 case 1: // conductor at specified potential
672 if (boundarytype2 == 0 || boundarytype2 == 4) {
673 // dielectric-conductor
674 InterfaceType[prim] = 1;
675 ApplPot[prim] = potential1;
676 } else if (boundarytype2 == 1) {
677 // conductor-conductor
678 if (fabs(potential1 - potential2) // same potential
679 < 1e-6 * (1 + fabs(potential1) + fabs(potential2))) {
680 printf("neBEMReadGeometry: identical potentials; skipped.\n");
681 printf("Primitive skipped: #%d\n", prim);
682 InterfaceType[prim] = 0;
683 } else {
684 // different potentials
685 printf("neBEMReadGeometry: different potentials; rejected.\n");
686 return -1;
687 }
688 } else {
689 // conductor-unknown
690 printf(
691 "neBEMReadGeometry: conductor at given potential; rejected.\n");
692 return -1;
693 }
694 break;
695
696 case 2: // conductor with a specified charge
697 if ((boundarytype2 == 0) || (boundarytype2 == 4)) {
698 // conductor-dielectric
699 InterfaceType[prim] = 2;
700 ApplCh[prim] = charge1;
701 } else {
702 printf("neBEMReadGeometry: charged conductor; rejected.\n");
703 return -1;
704 }
705 break;
706
707 case 3: // floating conductor (zero charge, perpendicular E)
708 if ((boundarytype2 == 0) || (boundarytype2 == 4)) {
709 // conductor-dielectric
710 InterfaceType[prim] = 3;
711 if (!NbFloatingConductors) // assuming only one floating conductor
712 NbFloatingConductors = 1; // in the system
713 } else {
714 printf("neBEMReadGeometry: floating conductor; rejected.\n");
715 return -1;
716 }
717 break;
718
719 case 4: // dielectric interface (plastic-plastic) without "manual" charge
720 if (boundarytype2 == 0) {
721 // dielectric-vacuum
722 // epsilon1 is self dielectric-constant
723 // epsilon2 is towards positive normal
724 InterfaceType[prim] = 4;
725 Lambda[prim] = (epsilon1 - epsilon2) / (epsilon1 + epsilon2);
726 // consistent with Bardhan's eqn 16 where (1 / (2*Lambda)) is used
727 } else if (boundarytype2 == 1) {
728 // dielectric-conductor
729 InterfaceType[prim] = 1; // conductor at known potential
730 ApplPot[prim] = potential2;
731 } else if (boundarytype2 == 2) {
732 // dielectric-conductor
733 InterfaceType[prim] = 2; // conductor with known charge
734 ApplCh[prim] = charge2;
735 } else if (boundarytype2 == 3) {
736 // dielectric-conductor
737 InterfaceType[prim] = 3; // conductor at floating potential
738 } else if (boundarytype2 == 4) {
739 // dielectric-dielectric
740 if (fabs(epsilon1 - epsilon2) < 1e-6 * (1 + fabs(epsilon1) + fabs(epsilon2))) {
741 // identical dielectrica
742 printf(
743 "neBEMReadGeometry: between identical dielectrica; skipd.\n");
744 printf("Primitive skipped: #%d\n", prim);
745 InterfaceType[prim] = 0;
746 } else {
747 // distinctly different dielectrica
748 // epsilon1 is self dielectric-constant
749 // epsilon2 towards positive normal
750 InterfaceType[prim] = 4;
751 Lambda[prim] = (epsilon1 - epsilon2) / (epsilon1 + epsilon2);
752 // consistent with Bardhan's paper (1 / Lambda)
753 }
754 } else if (boundarytype2 == 5) {
755 // dielectric-dielectric with charge
756 if (fabs(epsilon1 - epsilon2) // identical dielectrica
757 < 1e-6 * (1 + fabs(epsilon1) + fabs(epsilon2))) {
758 printf(
759 "neBEMReadGeometry: between identical dielectrica; skipped.\n");
760 printf("Primitive skipped: #%d\n", prim);
761 InterfaceType[prim] = 0;
762 } else {
763 // distinctly different dielectrica
764 InterfaceType[prim] = 5; // epsilon2 towards positive normal
765 ApplCh[prim] = charge2;
766 Lambda[prim] = (epsilon1 - epsilon2) / (epsilon1 + epsilon2);
767 }
768 } // if-else if boundarytypes 0 and 4
769 else { // dielectric-unknown
770 printf("neBEMReadGeometry: unknown dielectric; rejected.\n");
771 return -1;
772 }
773 break;
774
775 case 5: // dielectric with surface charge (plastic-gas, typically)
776 if (boundarytype2 == 0) { // dielectric-vacuum
777 InterfaceType[prim] = 5; // epsilon2 is towards +ve normal
778 ApplCh[prim] = charge1;
779 Lambda[prim] = (epsilon1 - epsilon2) / (epsilon1 + epsilon2);
780 // consistent with Bardhan's paper (1 / Lambda)
781 } else if (boundarytype2 == 4) { // dielectric-dielectric
782 if (fabs(epsilon1 - epsilon2) < 1e-6 * (1 + fabs(epsilon1) + fabs(epsilon2))) {
783 // identical dielectrica
784 printf(
785 "neBEMReadGeometry: between identical dielectrica; skipd.\n");
786 printf("Primitive skipped: #%d\n", prim);
787 InterfaceType[prim] = 0;
788 } else {
789 // distinctly different dielectrica
790 InterfaceType[prim] = 5; // epsilon2 towards positive normal
791 ApplCh[prim] = charge1;
792 Lambda[prim] = (epsilon1 - epsilon2) / (epsilon1 + epsilon2);
793 // consistent with Bardhan's paper (1 / Lambda)
794 }
795 } // if-else if boundarytypes 0 and 4
796 else {
797 printf(
798 "neBEMReadGeometry: charged dielectric adjacent to a conductor; "
799 "rejected.\n");
800 return -1;
801 }
802 break;
803
804 case 6: // symmetry boundary, E parallel
805 if (boundarytype2 == 0) {
806 InterfaceType[prim] = 6;
807 } else {
808 printf("neBEMReadGeometry: E-parallel symmetry; rejected.\n");
809 return -1;
810 }
811 break;
812
813 case 7: // symmetry boundary, E perpendicular
814 if (boundarytype2 == 0) {
815 InterfaceType[prim] = 7;
816 } else {
817 printf("neBEMReadGeometry: E-perpendicular symmetry; rejected.\n");
818 return -1;
819 }
820 break;
821
822 default:
823 printf("neBEMReadGeometry: Boundary type 1: %d\n", boundarytype1);
824 printf("neBEMReadGeometry: Boundary type 2: %d\n", boundarytype2);
825 printf("neBEMReadGeometry: out of range ... exiting.\n");
826 return -1;
827 } // switch boundarytype1 ends
828
830 printf(
831 "\tType: %d, ApplPot: %lg, Epsilon1: %lg, Epsilon2: %lg, Lambda: "
832 "%lg, ApplCh: %lg\n",
833 InterfaceType[prim], ApplPot[prim], Epsilon1[prim], Epsilon2[prim],
834 Lambda[prim], ApplCh[prim]);
835 }
836
837 // Read the periodicities
838 // Note that mirror has been taken care of separately below
839 // ix: PeriodicTypeX (1 - simple, 2 - mirror, 3 - axial, 4 - rotation)
840 // jx: PeriodicInX (Number of copies internal to neBEM)
841 // sx: XPeriod
842 // NOTE: A change in this part of the algorithm is likely to affect
843 // src/Solve/neBEM.c (LHMatrix) and
844 // src/Solve/ComputeProperties.c (PFAtPoint and WtPFAtPoint)
845 {
846 int ix, iy, iz;
847 int jx, jy, jz;
848 double sx, sy, sz;
849 fstatus = neBEMGetPeriodicities(prim, &ix, &jx, &sx, &iy, &jy, &sy, &iz,
850 &jz, &sz);
851 if (fstatus != 0) {
852 neBEMMessage("neBEMReadGeometry - neBEMGetPeriodicities");
853 return -1;
854 }
855 if (jx < 0) jx = 0;
856 if (jy < 0) jy = 0;
857 if (jz < 0) jz = 0;
858
859 PeriodicTypeX[prim] = ix;
860 PeriodicTypeY[prim] = iy;
861 PeriodicTypeZ[prim] = iz;
862 if (0) {
863 printf("For primitive: %d\n", prim);
864 printf("\tPeriodicTypeX: %d, PeriodicTypeY: %d, PeriodicTypeZ: %d\n",
865 ix, iy, iz);
866 printf("\tPeriodicInX: %d, PeriodicInY: %d, PeriodicInZ: %d\n", jx, jy,
867 jz);
868 printf("\tXPeriod: %lg, YPeriod: %lg, ZPeriod: %lg\n", sx, sy, sz);
869 }
870 if (ix > 0) {
871 // These checks need to be done separately. Otherwise, there is
872 // a possibility that non-zero values of PeriodicIn* and *Period
873 // are used throughout the code despite PeriodicType* is 0
874 PeriodicInX[prim] = jx;
875 XPeriod[prim] = sx;
876 } else {
877 PeriodicInX[prim] = 0;
878 XPeriod[prim] = 0.0;
879 }
880 if (iy > 0) {
881 PeriodicInY[prim] = jy;
882 YPeriod[prim] = sy;
883 } else {
884 PeriodicInY[prim] = 0;
885 YPeriod[prim] = 0.0;
886 }
887 if (iz > 0) {
888 PeriodicInZ[prim] = jz;
889 ZPeriod[prim] = sz;
890 } else {
891 PeriodicInZ[prim] = 0;
892 ZPeriod[prim] = 0.0;
893 }
894 } // read periodicity information
895
896 // Read mirror information
897 // Mirror can be in perpendicular to any of the three cartesian axes
898 // There can be more than one mirror at the same time
899 // MirrorType (1 - charge density -ve of original, equivalent to method of
900 // images, 2 - charge density same as original)
901 {
902 int ix, iy, iz;
903 int jx, jy, jz; // not used at present
904 double sx, sy, sz;
905 fstatus =
906 neBEMGetMirror(prim, &ix, &jx, &sx, &iy, &jy, &sy, &iz, &jz, &sz);
907 if (fstatus != 0) {
908 neBEMMessage("neBEMReadGeometry - neBEMGetMirror");
909 return -1;
910 }
911 if (jx < 0) jx = 0;
912 if (jy < 0) jy = 0;
913 if (jz < 0) jz = 0;
914
915 MirrorTypeX[prim] = ix;
916 MirrorTypeY[prim] = iy;
917 MirrorTypeZ[prim] = iz;
918 if (0) {
919 printf("For primitive: %d\n", prim);
920 printf("\tMirrorTypeX: %d, MirrorTypeY: %d, MirrorTypeZ: %d\n", ix, iy,
921 iz);
922 printf("\tNOT USED ==> MirrorInX: %d, MirrorInY: %d, MirrorInZ: %d\n",
923 jx, jy, jz);
924 printf("\tMirrorDistX: %lg, MirrorDistY: %lg, MirrorDistZ: %lg\n", sx,
925 sy, sz);
926 getchar();
927 }
928 if (ix > 0) {
929 // printf("neBEMReadGeometry: Mirror have been requested.\n");
930 MirrorDistXFromOrigin[prim] = sx; // assumed to pass through the origin
931 } else {
932 MirrorDistXFromOrigin[prim] = 0.0; // pass through the origin
933 }
934 if (iy > 0) {
935 // printf("neBEMReadGeometry: Mirror have been requested.\n");
936 MirrorDistYFromOrigin[prim] = sy;
937 } else {
938 MirrorDistYFromOrigin[prim] = 0.0;
939 }
940 if (iz > 0) {
941 // printf("neBEMReadGeometry: Mirror have been requested.\n");
942 MirrorDistZFromOrigin[prim] = sz;
943 } else {
944 MirrorDistZFromOrigin[prim] = 0.0;
945 }
946 } // read mirror information
947
948 // Information on bounding planes
949 // ixmin=0: lower x-plane does not exist
950 // ixmin=1: lower x-plane does exist
951 // cxmin: coordinate of lower x-plane
952 // vxmin: potential of lower x-plane
953 // Similar for ixmax, iymin, iymax, izmin, izmax
954 int ixmin, ixmax, iymin, iymax, izmin, izmax;
955 double cxmin, cxmax, cymin, cymax, czmin, czmax;
956 double vxmin, vxmax, vymin, vymax, vzmin, vzmax;
957 fstatus = neBEMGetBoundingPlanes(
958 &ixmin, &cxmin, &vxmin, &ixmax, &cxmax, &vxmax, &iymin, &cymin, &vymin,
959 &iymax, &cymax, &vymax, &izmin, &czmin, &vzmin, &izmax, &czmax, &vzmax);
960 if (fstatus != 0) {
961 neBEMMessage("neBEMReadGeometry - neBEMGetBoundingPlanes");
962 return -1;
963 }
964 BndPlaneInXMin[prim] = ixmin;
965 BndPlaneInXMax[prim] = ixmax;
966 BndPlaneInYMin[prim] = iymin;
967 BndPlaneInYMax[prim] = iymax;
968 BndPlaneInZMin[prim] = izmin;
969 BndPlaneInZMax[prim] = izmax;
970 if (ixmin) {
971 XBndPlaneInXMin[prim] = cxmin;
972 VBndPlaneInXMin[prim] = vxmin;
973 } else {
974 XBndPlaneInXMin[prim] = 0.0;
975 VBndPlaneInXMin[prim] = 0.0;
976 }
977 if (ixmax > 0) {
978 XBndPlaneInXMax[prim] = cxmax;
979 VBndPlaneInXMax[prim] = vxmax;
980 } else {
981 XBndPlaneInXMax[prim] = 0.0;
982 VBndPlaneInXMax[prim] = 0.0;
983 }
984 if (iymin > 0) {
985 YBndPlaneInYMin[prim] = cymin;
986 VBndPlaneInYMin[prim] = vymin;
987 } else {
988 YBndPlaneInYMin[prim] = 0.0;
989 VBndPlaneInYMin[prim] = 0.0;
990 }
991 if (iymax > 0) {
992 YBndPlaneInYMax[prim] = cymax;
993 VBndPlaneInYMax[prim] = vymax;
994 } else {
995 YBndPlaneInYMax[prim] = 0.0;
996 VBndPlaneInYMax[prim] = 0.0;
997 }
998 if (izmin > 0) {
999 ZBndPlaneInZMin[prim] = czmin;
1000 VBndPlaneInZMin[prim] = vzmin;
1001 } else {
1002 ZBndPlaneInZMin[prim] = 0.0;
1003 VBndPlaneInZMin[prim] = 0.0;
1004 }
1005 if (izmax > 0) {
1006 ZBndPlaneInZMax[prim] = czmax;
1007 VBndPlaneInZMax[prim] = vzmax;
1008 } else {
1009 ZBndPlaneInZMax[prim] = 0.0;
1010 VBndPlaneInZMax[prim] = 0.0;
1011 }
1012 } // loop over the primitives - major loop
1013
1014 VolMax = volmax; // Maximum nb of volumes in the problem
1015 volRef = ivector(0, VolMax); // variables to store volume related infomration
1016 volShape = ivector(0, VolMax);
1020 volCharge = dvector(0, VolMax);
1022 for (int volref = 0; volref <= VolMax; ++volref) {
1023 neBEMVolumeDescription(volref, &volShape[volref], &volMaterial[volref],
1024 &volEpsilon[volref], &volPotential[volref],
1025 &volCharge[volref], &volBoundaryType[volref]);
1026 if (dbgFn) {
1027 printf("volref: %d\n", volref);
1028 printf("shape: %d, material: %d\n", volShape[volref],
1029 volMaterial[volref]);
1030 printf("eps: %lg, pot: %lg\n", volEpsilon[volref], volPotential[volref]);
1031 printf("q: %lg, type: %d\n", volCharge[volref], volBoundaryType[volref]);
1032 }
1033 }
1034
1035 // Ignore unnecessary primitives from the final count
1036 // Ideally, all the removal conditions for a primitive should be checked in
1037 // one loop and the list should be updated in one single go.
1038 {
1039 for (int prim = 1; prim <= NbPrimitives; ++prim) {
1040 OrgnlToEffPrim[prim][0] = prim;
1041 OrgnlToEffPrim[prim][1] = prim;
1042 OrgnlToEffPrim[prim][2] = prim;
1043 }
1044
1045 { // Skip primitive
1046 // Remove skipped primitives having InterfaceType == 0.
1047 // Also remove primitives having too small dimensions.
1048 int NbSkipped = 0, effprim;
1049 double DVertex[4], minDVertex = 0.0; // maximum number of vertices is 4
1050 for (int prim = 1; prim <= NbPrimitives; ++prim) {
1051 effprim = prim - NbSkipped;
1052
1053 // Check dimensions of the primitive
1054 for (int vert = 0; vert < NbVertices[prim] - 1; ++vert) {
1055 DVertex[vert] =
1056 sqrt(((XVertex[prim][vert + 1] - XVertex[prim][vert]) *
1057 (XVertex[prim][vert + 1] - XVertex[prim][vert])) +
1058 ((YVertex[prim][vert + 1] - YVertex[prim][vert]) *
1059 (YVertex[prim][vert + 1] - YVertex[prim][vert])) +
1060 ((ZVertex[prim][vert + 1] - ZVertex[prim][vert]) *
1061 (ZVertex[prim][vert + 1] - ZVertex[prim][vert])));
1062 if (vert == 0)
1063 minDVertex = DVertex[vert];
1064 else {
1065 if (DVertex[vert] < minDVertex) minDVertex = DVertex[vert];
1066 }
1067 }
1068
1069 if ((InterfaceType[prim]) && (minDVertex > MINDIST)) {
1070 OrgnlToEffPrim[prim][1] = effprim;
1071 OrgnlToEffPrim[prim][2] = effprim;
1072 PrimType[effprim] = PrimType[prim];
1073 NbVertices[effprim] = NbVertices[prim];
1074 for (int vert = 0; vert < NbVertices[effprim]; ++vert) {
1075 XVertex[effprim][vert] = XVertex[prim][vert];
1076 YVertex[effprim][vert] = YVertex[prim][vert];
1077 ZVertex[effprim][vert] = ZVertex[prim][vert];
1078 }
1079 if (PrimType[effprim] == 2) // wire
1080 {
1081 XNorm[effprim] = 0.0; // modulus not 1 - an absurd trio!
1082 YNorm[effprim] = 0.0;
1083 ZNorm[effprim] = 0.0;
1084 Radius[effprim] = Radius[prim];
1085 }
1086 if ((PrimType[effprim] == 3) || (PrimType[effprim] == 4)) {
1087 XNorm[effprim] = XNorm[prim];
1088 YNorm[effprim] = YNorm[prim];
1089 ZNorm[effprim] = ZNorm[prim];
1090 Radius[effprim] = 0.0; // absurd radius!
1091 }
1092 VolRef1[effprim] = VolRef1[prim];
1093 VolRef2[effprim] = VolRef2[prim];
1094
1095 InterfaceType[effprim] = InterfaceType[prim];
1096 Epsilon1[effprim] = Epsilon1[prim];
1097 Epsilon2[effprim] = Epsilon2[prim];
1098 Lambda[effprim] = Lambda[prim];
1099 ApplPot[effprim] = ApplPot[prim];
1100 ApplCh[effprim] = ApplCh[prim];
1101 PeriodicTypeX[effprim] = PeriodicTypeX[prim];
1102 PeriodicTypeY[effprim] = PeriodicTypeY[prim];
1103 PeriodicTypeZ[effprim] = PeriodicTypeZ[prim];
1104 PeriodicInX[effprim] = PeriodicInX[prim];
1105 PeriodicInY[effprim] = PeriodicInY[prim];
1106 PeriodicInZ[effprim] = PeriodicInZ[prim];
1107 XPeriod[effprim] = XPeriod[prim];
1108 YPeriod[effprim] = YPeriod[prim];
1109 ZPeriod[effprim] = ZPeriod[prim];
1110 MirrorTypeX[effprim] = MirrorTypeX[prim];
1111 MirrorTypeY[effprim] = MirrorTypeY[prim];
1112 MirrorTypeZ[effprim] = MirrorTypeZ[prim];
1116 BndPlaneInXMin[effprim] = BndPlaneInXMin[prim];
1117 BndPlaneInXMax[effprim] = BndPlaneInXMax[prim];
1118 BndPlaneInYMin[effprim] = BndPlaneInYMin[prim];
1119 BndPlaneInYMax[effprim] = BndPlaneInYMax[prim];
1120 BndPlaneInZMin[effprim] = BndPlaneInZMin[prim];
1121 BndPlaneInZMax[effprim] = BndPlaneInZMax[prim];
1122 XBndPlaneInXMin[effprim] = XBndPlaneInXMin[prim];
1123 XBndPlaneInXMax[effprim] = XBndPlaneInXMax[prim];
1124 YBndPlaneInYMin[effprim] = YBndPlaneInYMin[prim];
1125 YBndPlaneInYMax[effprim] = YBndPlaneInYMax[prim];
1126 ZBndPlaneInZMin[effprim] = ZBndPlaneInZMin[prim];
1127 ZBndPlaneInZMax[effprim] = ZBndPlaneInZMax[prim];
1128 VBndPlaneInXMin[effprim] = VBndPlaneInXMin[prim];
1129 VBndPlaneInXMax[effprim] = VBndPlaneInXMax[prim];
1130 VBndPlaneInYMin[effprim] = VBndPlaneInYMin[prim];
1131 VBndPlaneInYMax[effprim] = VBndPlaneInYMax[prim];
1132 VBndPlaneInZMin[effprim] = VBndPlaneInZMin[prim];
1133 VBndPlaneInZMax[effprim] = VBndPlaneInZMax[prim];
1134 } // InterfaceType
1135 else {
1136 OrgnlToEffPrim[prim][1] = 0; // removed from the list
1137 OrgnlToEffPrim[prim][2] = 0;
1138 ++NbSkipped;
1139 if (DebugLevel == 101) {
1140 printf("Skipped primitive %d, InterfaceType: %d, minDVertex: %lg\n",
1141 prim, InterfaceType[prim], minDVertex);
1142 }
1143 }
1144 } // loop over primitives to remove the skipped primitives
1145 NbPrimitives -= NbSkipped;
1146 printf("Number of primitives skipped: %d, Effective NbPrimitives: %d\n",
1147 NbSkipped, NbPrimitives);
1148 } // Skip primitives
1149
1150 /*
1151 { // look for overlaps among primitives
1152 // Also remove primitives so that periodicities do not lead to overlapped
1153 // (full or partial?) primitives.
1154 // Check overlap of primitives due to repetition
1155 // Algorithm:
1156 // 1. Consider those primitives which have normal along X.
1157 // 2. Consider their initial positions in X.
1158 // 3. Find out their positions after one positive repetition in both X and
1159 Y.
1160 // and check which of them fall on already existing primitives having normal
1161 // along X.
1162 // A. Similarly for one negative repetition.
1163 // II. Same check for all the primitives having normal along Y will have to
1164 be
1165 // carried out.
1166 // Similar approach could be useful for an arbitrary primitive if we work in
1167 // the local coordinate system of the primitive and check all other
1168 primitives
1169 // in that LCS.
1170 // Complications are likely to arise due to corner and edge overlaps which
1171 will
1172 // happen to all neighbouring primitives on the same plane.
1173 int Overlap[NbPrimitives+1][NbPrimitives+1]; // this C++ syntax seems to
1174 work!
1175 // int **Overlap;
1176 // Overlap = imatrix(1, NbPrimitives, 1, NbPrimitives);
1177 for(unsigned int prim = 1; prim <= NbPrimitives; ++prim)
1178 for(unsigned int chkprim = 1;
1179 chkprim <= NbPrimitives; ++chkprim)
1180 Overlap[prim][chkprim] = 0;
1181
1182 for(unsigned int prim = 1; prim <= NbPrimitives; ++prim)
1183 {
1184 if(dbgFn)
1185 printf("\nNew XNorm[%d]: %lg: ", prim, XNorm[prim]);
1186
1187 if(fabs(fabs(XNorm[prim]) - 1.0) >= 1.0e-12) // primitive || to
1188 YZ plane continue;
1189
1190 for(unsigned int chkprim = prim+1;
1191 chkprim <= NbPrimitives; ++chkprim)
1192 {
1193 if(dbgFn)
1194 printf("XNorm[%d]: %lg, ", chkprim, XNorm[chkprim]);
1195
1196 if(fabs(fabs(XNorm[chkprim]) - 1.0) >= 1.0e-12) //
1197 primitive || to YZ plane continue;
1198
1199 if(fabs(XVertex[prim][0] - XVertex[chkprim][0]) <= 1.0e-12)
1200 { // same plane; check YZ vertices for full /
1201 partial overlap double smallY = YVertex[prim][0]; double bigY =
1202 YVertex[prim][0]; double smallZ = ZVertex[prim][0]; double bigZ =
1203 ZVertex[prim][0]; for(unsigned int vert = 1; vert<=
1204 NbVertices[prim]; ++vert)
1205 {
1206 if(smallY > YVertex[prim][vert])
1207 smallY = YVertex[prim][vert];
1208 if(bigY < YVertex[prim][vert])
1209 bigY = YVertex[prim][vert];
1210 if(smallZ > ZVertex[prim][vert])
1211 smallZ = ZVertex[prim][vert];
1212 if(bigZ < ZVertex[prim][vert])
1213 bigZ = ZVertex[prim][vert];
1214 }
1215
1216 double smallcY = YVertex[chkprim][0];
1217 double bigcY = YVertex[chkprim][0];
1218 double smallcZ = ZVertex[chkprim][0];
1219 double bigcZ = ZVertex[chkprim][0];
1220 for(unsigned int vert = 1;
1221 vert<= NbVertices[chkprim]; ++vert)
1222 {
1223 if(smallcY > YVertex[chkprim][vert])
1224 smallcY = YVertex[chkprim][vert];
1225 if(bigcY < YVertex[chkprim][vert])
1226 bigcY = YVertex[chkprim][vert];
1227 if(smallcZ > ZVertex[chkprim][vert])
1228 smallcZ = ZVertex[chkprim][vert];
1229 if(bigcZ < ZVertex[chkprim][vert])
1230 bigcZ = ZVertex[chkprim][vert];
1231 }
1232
1233 // Because of the rearrangement, the big and small
1234 vertices are now
1235 // distinct and unambiguous. Comparisons can now be
1236 made a lot more
1237 // easily.
1238 // One important assumption is implicit,however.
1239 // The edges of the primitive are assumed to
1240 parallel to the
1241 // Y and Z axes. We are assuming rectangular
1242 primitives as a
1243 // consequence.
1244 // One to one match of all the vertices
1245 // Subsequent checks should be carried out only if
1246 Overlap is still 0
1247
1248 // No overlap
1249 if(((smallcY > smallY) && (smallcY > bigY))
1250 && ((bigcY > smallY) && (bigcY >
1251 bigY))) Overlap[prim][chkprim] = -1; // no Y overlap
1252 if(Overlap[prim][chkprim] != 0) continue;
1253 if(((smallcZ > smallZ) && (smallcZ > bigZ))
1254 && ((bigcZ > smallZ) && (bigcZ >
1255 bigZ))) Overlap[prim][chkprim] = -2; // no Z overlap
1256 if(Overlap[prim][chkprim] != 0) continue;
1257
1258 // Check for shared corners
1259 if((fabs(bigcY-smallY) <= MINDIST)
1260 && (fabs(bigcZ-smallZ) <= MINDIST))
1261 Overlap[prim][chkprim] = -3;
1262 if(Overlap[prim][chkprim] != 0) continue;
1263 if((fabs(bigcY-smallY) <= MINDIST)
1264 && (fabs(smallcZ-bigZ) <= MINDIST))
1265 Overlap[prim][chkprim] = -4;
1266 if(Overlap[prim][chkprim] != 0) continue;
1267 if((fabs(smallcY-bigY) <= MINDIST)
1268 && (fabs(smallcZ-bigZ) <= MINDIST))
1269 Overlap[prim][chkprim] = -5;
1270 if(Overlap[prim][chkprim] != 0) continue;
1271 if((fabs(smallcY-bigY) <= MINDIST)
1272 && (fabs(bigcZ-bigZ) <= MINDIST))
1273 Overlap[prim][chkprim] = -6;
1274 if(Overlap[prim][chkprim] != 0) continue;
1275
1276 // Check for completely / partially shared edges
1277 if((fabs(bigcY-smallY) <= MINDIST)
1278 && (((smallcZ > smallZ) && (smallcZ
1279 < bigZ)) // = implies corner
1280 || ((bigcZ > smallZ) && (bigcZ <
1281 bigZ)))) // = implies corner Overlap[prim][chkprim] = -7;
1282 if(Overlap[prim][chkprim] != 0) continue;
1283 if((fabs(smallcZ-bigZ) <= MINDIST)
1284 && (((smallcY > smallY) && (smallcY
1285 < bigY)) // = implies corner
1286 || ((bigcY > smallY) && (bigcY <
1287 bigY)))) // = implies corner Overlap[prim][chkprim] = -8;
1288 if(Overlap[prim][chkprim] != 0) continue;
1289 if((fabs(smallcY-bigY) <= MINDIST)
1290 && (((smallcZ > smallZ) && (smallcZ
1291 < bigZ)) // = implies corner
1292 || ((bigcZ > smallZ) && (bigcZ <
1293 bigZ)))) // = implies corner Overlap[prim][chkprim] = -9;
1294 if(Overlap[prim][chkprim] != 0) continue;
1295 if((fabs(bigcZ-smallZ) <= MINDIST)
1296 && (((smallcY > smallY) && (smallcY
1297 < bigY)) // = implies corner
1298 || ((bigcY > smallY) && (bigcY <
1299 bigY)))) // = implies corner Overlap[prim][chkprim] = -10;
1300 if(Overlap[prim][chkprim] != 0) continue;
1301
1302 // One-to-one overlap - case 0
1303 if((fabs(smallcY-smallY) <= MINDIST)
1304 && (fabs(bigcY-bigY) <= MINDIST)
1305 && (fabs(smallcZ-smallZ) <= MINDIST)
1306 && (fabs(bigcZ-bigZ) <= MINDIST))
1307 { // chkprim has a one-to-one overlap with
1308 prim Overlap[prim][chkprim] = 1;
1309 }
1310 if(Overlap[prim][chkprim] != 0) continue;
1311
1312 // checkprim is completely contained within prim -
1313 case 1 if((smallcY > smallY) && (smallcY < bigY)
1314 && (smallcZ > smallZ) && (smallcZ <
1315 bigZ)
1316 && (bigcY > smallY) && (bigcY <
1317 bigY)
1318 && (bigcZ > smallZ) && (bigcZ <
1319 bigZ)) { // chkprim completely contained within prim
1320 Overlap[prim][chkprim] = 2;
1321 }
1322 if(Overlap[prim][chkprim] != 0) continue;
1323
1324 // checkprim completely contains prim - case 2
1325 if((smallY > smallcY) && (smallY < bigcY)
1326 && (smallZ > smallcZ) && (smallZ <
1327 bigcZ)
1328 && (bigY > smallcY) && (bigY <
1329 bigcY)
1330 && (bigZ > smallcZ) && (bigZ <
1331 bigcZ)) { // chkprim completely contains prim Overlap[prim][chkprim] = 3;
1332 }
1333 if(Overlap[prim][chkprim] != 0) continue;
1334
1335 // checkprim completely contained within prim in Z
1336 direction - case 3 if((smallcZ > smallZ) && (bigcZ < bigZ))
1337 {
1338 if((smallcY > smallY) && (smallcY < bigY)
1339 && (bigcY > smallY) &&
1340 (bigcY > bigY)) { // but not Y: bigY outside prim Overlap[prim][chkprim] =
1341 4;
1342 }
1343 if(Overlap[prim][chkprim] != 0) continue;
1344 if((smallcY < smallY) && (smallcY < bigY)
1345 && (bigcY > smallY) &&
1346 (bigcY < bigY)) { // but not Y: smallY outside prim Overlap[prim][chkprim]
1347 = 5;
1348 }
1349 if(Overlap[prim][chkprim] != 0) continue;
1350 if((smallcY < smallY) && (smallcY < bigY)
1351 && (bigcY > smallY) &&
1352 (bigcY > bigY)) { // but not Y: chkprim wider than prim
1353 Overlap[prim][chkprim] = 6;
1354 }
1355 if(Overlap[prim][chkprim] != 0) continue;
1356 }
1357
1358 // checkprim completely contained within prim in Y
1359 direction - case 4 if((smallcY > smallY) && (bigcY < bigY))
1360 {
1361 if((smallcZ > smallZ) && (smallcZ < bigZ)
1362 // bigcZ outside
1363 && (bigcZ > smallZ) &&
1364 (bigcZ > bigZ)) { // but not Z: bigZ outside prim Overlap[prim][chkprim] =
1365 7;
1366 }
1367 if(Overlap[prim][chkprim] != 0) continue;
1368 if((smallcZ < smallZ) && (smallcZ < bigZ)
1369 // smallcZ outside
1370 && (bigcZ > smallZ) &&
1371 (bigcZ < bigZ)) { // but not Z: smallZ outside prim Overlap[prim][chkprim]
1372 = 8;
1373 }
1374 if(Overlap[prim][chkprim] != 0) continue;
1375 if((smallcZ < smallZ) && (smallcZ < bigZ)
1376 // smallcZ outside
1377 && (bigcZ > smallZ) &&
1378 (bigcZ > bigZ)) // bigcZ outside { // but not Z: chkprim taller than prim
1379 Overlap[prim][chkprim] = 9;
1380 }
1381 if(Overlap[prim][chkprim] != 0) continue;
1382 }
1383
1384 // checkprim completely contains prim in Z direction
1385 - case 5 if((smallcZ < smallZ) && (bigcZ > bigZ))
1386 {
1387 if((smallY > smallcY) && (smallY < bigcY)
1388 // bigY outside
1389 && (bigY > smallcY) && (bigY
1390 > bigcY)) { // but not Y: bigY outside chkprim Overlap[prim][chkprim] = 10;
1391 }
1392 if(Overlap[prim][chkprim] != 0) continue;
1393 if((smallY < smallcY) && (smallY < bigcY)
1394 // smallY outside
1395 && (bigY > smallcY) && (bigY
1396 < bigcY)) { // but not Y: smallY outside chkprim Overlap[prim][chkprim] =
1397 11;
1398 }
1399 if(Overlap[prim][chkprim] != 0) continue;
1400 if((smallY < smallcY) && (bigY > bigcY))
1401 // outside chkprim extent { // but not Y: prim wider than chkprim
1402 Overlap[prim][chkprim] = 12;
1403 }
1404 if(Overlap[prim][chkprim] != 0) continue;
1405 }
1406
1407 // checkprim completely contains prim in Y direction
1408 - case 6 if((smallcY < smallY) && (bigcY > bigY))
1409 {
1410 if((smallZ > smallcZ) && (smallZ < bigcZ)
1411 // bigZ outside
1412 && (bigZ > smallcZ) && (bigZ
1413 > bigcZ)) { // but not Z: bigZ outside chkprim Overlap[prim][chkprim] = 13;
1414 }
1415 if(Overlap[prim][chkprim] != 0) continue;
1416 if((smallZ < smallcZ) && (smallZ < bigcZ)
1417 // smallZ outside
1418 && (bigZ > smallcZ) && (bigZ
1419 < bigcZ)) { // but not Z: smallZ outside chkprim Overlap[prim][chkprim] =
1420 14;
1421 }
1422 if(Overlap[prim][chkprim] != 0) continue;
1423 if((smallZ < smallcZ) && (smallZ < bigcZ)
1424 // smallZ and bigZ outside
1425 && (bigZ > smallcZ) && (bigZ
1426 < bigcZ)) { // but not Z: prim taller than chkprim Overlap[prim][chkprim] =
1427 14;
1428 }
1429 if(Overlap[prim][chkprim] != 0) continue;
1430 }
1431
1432 // check cases where only one corner of chkprim is
1433 within prim
1434 // all the other three being beyond.
1435 if((smallcY > smallY) && (smallcY < bigY)
1436 && (smallcZ > smallZ) && (smallcZ <
1437 bigZ)
1438 && (bigcY > bigY) && (bigcZ > bigZ))
1439 { // only smallcY,smallcZ corner within
1440 prim Overlap[prim][chkprim] = 15;
1441 }
1442 if(Overlap[prim][chkprim] != 0) continue;
1443 if((bigcY > smallY) && (bigcY < bigY)
1444 && (smallcZ > smallZ) && (smallcZ <
1445 bigZ)
1446 && (smallcY < smallY) && (bigcZ >
1447 bigZ)) { // only bigcY,smallcZ corner within prim Overlap[prim][chkprim]
1448 = 16;
1449 }
1450 if(Overlap[prim][chkprim] != 0) continue;
1451 if((smallcY > smallY) && (smallcY < bigY)
1452 && (bigcZ > smallZ) && (bigcZ <
1453 bigZ)
1454 && (bigcY > bigY) && (smallcZ <
1455 smallZ)) { // only smallcY,bigcZ corner within prim Overlap[prim][chkprim]
1456 = 17;
1457 }
1458 if(Overlap[prim][chkprim] != 0) continue;
1459 if((bigcY > smallY) && (bigcY < bigY)
1460 && (bigcZ > smallZ) && (bigcZ <
1461 bigZ)
1462 && (smallcY < smallY) && (smallcZ <
1463 smallZ)) { // only bigcY,bigcZ corner within prim Overlap[prim][chkprim] =
1464 18;
1465 }
1466 if(Overlap[prim][chkprim] != 0) continue;
1467
1468 // check cases where only one corner of prim is
1469 within chkprim
1470 // all the other three being beyond.
1471 if((smallY > smallcY) && (smallY < bigcY)
1472 && (smallZ > smallcZ) && (smallZ <
1473 bigcZ)
1474 && (bigcY < bigY) && (bigcZ < bigZ))
1475 { // only smallY,smallZ corner within prim
1476 Overlap[prim][chkprim] = 19;
1477 }
1478 if(Overlap[prim][chkprim] != 0) continue;
1479 if((bigY > smallcY) && (bigY < bigcY)
1480 && (smallZ > smallcZ) && (smallZ <
1481 bigcZ)
1482 && (smallcY > smallY) && (bigcZ <
1483 bigZ)) { // only bigY,smallZ corner within prim Overlap[prim][chkprim] =
1484 20;
1485 }
1486 if(Overlap[prim][chkprim] != 0) continue;
1487 if((smallY > smallcY) && (smallY < bigcY)
1488 && (bigZ > smallcZ) && (bigZ <
1489 bigcZ)
1490 && (bigcY < bigY) && (smallcZ >
1491 smallZ)) { // only smallY,bigZ corner within prim Overlap[prim][chkprim] =
1492 21;
1493 }
1494 if(Overlap[prim][chkprim] != 0) continue;
1495 if((bigY > smallcY) && (bigY < bigcY)
1496 && (bigZ > smallcZ) && (bigZ <
1497 bigcZ)
1498 && (smallcY > smallY) && (smallcZ >
1499 smallZ)) { // only bigY,bigZ corner within prim Overlap[prim][chkprim] =
1500 22;
1501 }
1502 if(Overlap[prim][chkprim] != 0) continue;
1503 } // same plane; check overlap
1504
1505 } // loop over check primitives
1506 } // loop over primitives
1507
1508 if(dbgFn) printf("\n");
1509 for(unsigned int prim = 1; prim <= NbPrimitives; ++prim)
1510 for(unsigned int chkprim = prim+1;
1511 chkprim <= NbPrimitives; ++chkprim)
1512 {
1513 printf("prim: %d, XNorm: %lg, chkprim: %d, XNorm: %lg,
1514 Overlap: %d\n", prim, XNorm[prim], chkprim, XNorm[chkprim],
1515 Overlap[prim][chkprim]);
1516 }
1517 } // look for overlapped primitives
1518 */
1519
1520 { // remove primitives, as specified in a user supplied input file
1521 int NbRmPrims;
1522 FILE *rmprimFile = fopen("neBEMRmPrim.inp", "r");
1523 if (rmprimFile == NULL) {
1524 printf("neBEMRmPrim.inp absent ... assuming defaults ...\n");
1525 NbRmPrims = 0;
1526 } else {
1527 fscanf(rmprimFile, "NbRmPrims: %d\n", &NbRmPrims);
1528 if (NbRmPrims) {
1529 int tint;
1530#ifdef __cplusplus
1531 std::vector<double> rmXNorm(NbRmPrims + 1, 0.);
1532 std::vector<double> rmYNorm(NbRmPrims + 1, 0.);
1533 std::vector<double> rmZNorm(NbRmPrims + 1, 0.);
1534 std::vector<double> rmXVert(NbRmPrims + 1, 0.);
1535 std::vector<double> rmYVert(NbRmPrims + 1, 0.);
1536 std::vector<double> rmZVert(NbRmPrims + 1, 0.);
1537#else
1538 double rmXNorm[NbRmPrims + 1], rmYNorm[NbRmPrims + 1];
1539 double rmZNorm[NbRmPrims + 1];
1540 double rmXVert[NbRmPrims + 1], rmYVert[NbRmPrims + 1];
1541 double rmZVert[NbRmPrims + 1];
1542#endif
1543 for (int rmprim = 1; rmprim <= NbRmPrims; ++rmprim) {
1544 fscanf(rmprimFile, "Prim: %d\n", &tint);
1545 fscanf(rmprimFile, "rmXNorm: %le\n", &rmXNorm[rmprim]);
1546 fscanf(rmprimFile, "rmYNorm: %le\n", &rmYNorm[rmprim]);
1547 fscanf(rmprimFile, "rmZNorm: %le\n", &rmZNorm[rmprim]);
1548 fscanf(rmprimFile, "rmXVert: %le\n", &rmXVert[rmprim]);
1549 fscanf(rmprimFile, "rmYVert: %le\n", &rmYVert[rmprim]);
1550 fscanf(rmprimFile, "rmZVert: %le\n", &rmZVert[rmprim]);
1551 printf(
1552 "rmprim: %d, rmXNorm: %lg, rmYNorm: %lg, rmZNorm: %lg, "
1553 "rmXVert: %lg, rmYVert: %lg, rmZVert: %lg\n",
1554 rmprim, rmXNorm[rmprim], rmYNorm[rmprim], rmZNorm[rmprim],
1555 rmXVert[rmprim], rmYVert[rmprim], rmZVert[rmprim]);
1556 }
1557#ifdef __cplusplus
1558 std::vector<int> remove(NbPrimitives + 1, 0);
1559#else
1560 int remove[NbPrimitives + 1];
1561#endif
1562 // Check updated prim list
1563 for (int prim = 1; prim <= NbPrimitives; ++prim) {
1564 remove[prim] = 0;
1565 if (dbgFn) {
1566 printf("\n\nprim: %d, XVertex: %lg, YVertex: %lg, ZVertex: %lg\n",
1567 prim, XVertex[prim][0], YVertex[prim][0],
1568 ZVertex[prim][0]);
1569 printf("XNorm: %lg, YNorm: %lg, ZNorm: %lg\n", XNorm[prim],
1570 YNorm[prim], ZNorm[prim]);
1571 }
1572
1573 for (int rmprim = 1; rmprim <= NbRmPrims; ++rmprim) {
1574 if (dbgFn) {
1575 printf(
1576 "rmprim: %d, rmXVertex: %lg, rmYVertex: %lg, rmZVertex: "
1577 "%lg\n",
1578 rmprim, rmXVert[rmprim], rmYVert[rmprim], rmZVert[rmprim]);
1579 printf("rmXNorm: %lg, rmYNorm: %lg, rmZNorm: %lg\n",
1580 rmXNorm[rmprim], rmYNorm[rmprim], rmZNorm[rmprim]);
1581 }
1582
1583 // check the normal
1584 if ((fabs(fabs(XNorm[prim]) - fabs(rmXNorm[rmprim])) <=
1585 MINDIST) &&
1586 (fabs(fabs(YNorm[prim]) - fabs(rmYNorm[rmprim])) <=
1587 MINDIST) &&
1588 (fabs(fabs(ZNorm[prim]) - fabs(rmZNorm[rmprim])) <=
1589 MINDIST)) { // prim and rmprim are parallel
1590 // coplanarity check to be implemented later.
1591 // For the time-being, we will assume that the planes to be
1592 // removed have their normals parallel to a given axis. So, we
1593 // only check that and remove the primitive if the distace along
1594 // that axis match. Possible pitfall => the primitives may be
1595 // coplanar but non-overlapping!
1596 if (fabs(fabs(XNorm[prim]) - 1.0) <= 1.0e-12) {
1597 // primitive || to YZ
1598 if (fabs(XVertex[prim][0] - rmXVert[rmprim]) <= MINDIST) {
1599 remove[prim] = 1;
1600 }
1601 }
1602 if (fabs(fabs(YNorm[prim]) - 1.0) <= 1.0e-12) {
1603 // primitive || to XZ
1604 if (fabs(YVertex[prim][0] - rmYVert[rmprim]) <= MINDIST) {
1605 remove[prim] = 1;
1606 }
1607 }
1608 if (fabs(fabs(ZNorm[prim]) - 1.0) <= 1.0e-12) {
1609 // primitive || to XY
1610 if (fabs(ZVertex[prim][0] - rmZVert[rmprim]) <= MINDIST) {
1611 remove[prim] = 1;
1612 }
1613 }
1614 } // case where prim and rmprim are parallel
1615 if (dbgFn) {
1616 printf("prim: %d, rmprim: %d, remove: %d\n", prim, rmprim,
1617 remove[prim]);
1618 }
1619 if (remove[prim] == 1)
1620 break; // once removed, no point checking others
1621 } // for rmprim - loop over all removal specification
1622
1623 } // for prim loop over all primitives
1624
1625 int NbRemoved = 0;
1626 FILE *fprrm = fopen("RmPrims.info", "w");
1627 if (fprrm == NULL) {
1628 printf(
1629 "error opening RmPrims.info file in write mode ... "
1630 "returning\n");
1631 return (-1);
1632 }
1633 // Note that some of the original primitives have already been removed
1634 // based on interface and dimension considerations
1635 int orgnlNb = 0;
1636 for (int prim = 1; prim <= NbPrimitives; ++prim) {
1637 // identify primitive number in the original list
1638 for (int orgnlprim = 1; orgnlprim <= OrgnlNbPrimitives;
1639 ++orgnlprim) {
1640 if (OrgnlToEffPrim[orgnlprim][1] ==
1641 prim) // number updated for intrfc
1642 {
1643 orgnlNb = orgnlprim;
1644 break;
1645 }
1646 } // loop for finding out its position in the original list
1647
1648 if (remove[prim] == 1) {
1649 ++NbRemoved;
1650 OrgnlToEffPrim[orgnlNb][2] = 0;
1651 fprintf(fprrm, "NbRemoved: %d, Removed primitive: %d\n",
1652 NbRemoved, prim);
1653 fprintf(fprrm, "PrimType: %d\n", PrimType[prim]);
1654 fprintf(fprrm, "NbVertices: %d\n", NbVertices[prim]);
1655 for (int vert = 0; vert < NbVertices[prim]; ++vert) {
1656 fprintf(fprrm, "Vertx %d: %lg, %lg, %lg\n", vert,
1657 XVertex[prim][vert], YVertex[prim][vert],
1658 ZVertex[prim][vert]);
1659 }
1660 fprintf(fprrm, "Normals: %lg, %lg, %lg\n", XNorm[prim],
1661 YNorm[prim], ZNorm[prim]);
1662 continue;
1663 } // if remove
1664 else { // keep this one in the updated list of primitives
1665 int effprim = prim - NbRemoved;
1666
1667 OrgnlToEffPrim[orgnlNb][2] = effprim;
1668 PrimType[effprim] = PrimType[prim];
1669 NbVertices[effprim] = NbVertices[prim];
1670 for (int vert = 0; vert < NbVertices[effprim]; ++vert) {
1671 XVertex[effprim][vert] = XVertex[prim][vert];
1672 YVertex[effprim][vert] = YVertex[prim][vert];
1673 ZVertex[effprim][vert] = ZVertex[prim][vert];
1674 }
1675 if (PrimType[effprim] == 2) {
1676 // wire
1677 XNorm[effprim] = 0.0; // modulus not 1 - an absurd trio!
1678 YNorm[effprim] = 0.0;
1679 ZNorm[effprim] = 0.0;
1680 Radius[effprim] = Radius[prim];
1681 }
1682 if ((PrimType[effprim] == 3) || (PrimType[effprim] == 4)) {
1683 XNorm[effprim] = XNorm[prim];
1684 YNorm[effprim] = YNorm[prim];
1685 ZNorm[effprim] = ZNorm[prim];
1686 Radius[effprim] = 0.0; // absurd radius!
1687 }
1688 VolRef1[effprim] = VolRef1[prim];
1689 VolRef2[effprim] = VolRef2[prim];
1690
1691 InterfaceType[effprim] = InterfaceType[prim];
1692 Epsilon1[effprim] = Epsilon1[prim];
1693 Epsilon2[effprim] = Epsilon2[prim];
1694 Lambda[effprim] = Lambda[prim];
1695 ApplPot[effprim] = ApplPot[prim];
1696 ApplCh[effprim] = ApplCh[prim];
1697 PeriodicTypeX[effprim] = PeriodicTypeX[prim];
1698 PeriodicTypeY[effprim] = PeriodicTypeY[prim];
1699 PeriodicTypeZ[effprim] = PeriodicTypeZ[prim];
1700 PeriodicInX[effprim] = PeriodicInX[prim];
1701 PeriodicInY[effprim] = PeriodicInY[prim];
1702 PeriodicInZ[effprim] = PeriodicInZ[prim];
1703 XPeriod[effprim] = XPeriod[prim];
1704 YPeriod[effprim] = YPeriod[prim];
1705 ZPeriod[effprim] = ZPeriod[prim];
1706 MirrorTypeX[effprim] = MirrorTypeX[prim];
1707 MirrorTypeY[effprim] = MirrorTypeY[prim];
1708 MirrorTypeZ[effprim] = MirrorTypeZ[prim];
1712 BndPlaneInXMin[effprim] = BndPlaneInXMin[prim];
1713 BndPlaneInXMax[effprim] = BndPlaneInXMax[prim];
1714 BndPlaneInYMin[effprim] = BndPlaneInYMin[prim];
1715 BndPlaneInYMax[effprim] = BndPlaneInYMax[prim];
1716 BndPlaneInZMin[effprim] = BndPlaneInZMin[prim];
1717 BndPlaneInZMax[effprim] = BndPlaneInZMax[prim];
1718 XBndPlaneInXMin[effprim] = XBndPlaneInXMin[prim];
1719 XBndPlaneInXMax[effprim] = XBndPlaneInXMax[prim];
1720 YBndPlaneInYMin[effprim] = YBndPlaneInYMin[prim];
1721 YBndPlaneInYMax[effprim] = YBndPlaneInYMax[prim];
1722 ZBndPlaneInZMin[effprim] = ZBndPlaneInZMin[prim];
1723 ZBndPlaneInZMax[effprim] = ZBndPlaneInZMax[prim];
1724 VBndPlaneInXMin[effprim] = VBndPlaneInXMin[prim];
1725 VBndPlaneInXMax[effprim] = VBndPlaneInXMax[prim];
1726 VBndPlaneInYMin[effprim] = VBndPlaneInYMin[prim];
1727 VBndPlaneInYMax[effprim] = VBndPlaneInYMax[prim];
1728 VBndPlaneInZMin[effprim] = VBndPlaneInZMin[prim];
1729 VBndPlaneInZMax[effprim] = VBndPlaneInZMax[prim];
1730 } // else remove == 0
1731 } // loop over primitives to remove the primitives tagged to be
1732 // removed
1733 fclose(fprrm);
1734
1735 NbPrimitives -= NbRemoved;
1736 printf(
1737 "Number of primitives removed: %d, Effective NbPrimitives: %d\n",
1738 NbRemoved, NbPrimitives);
1739 fflush(stdout);
1740 } // if NbRmPrimes true, implying primitives need to be removed
1741 fclose(rmprimFile);
1742 } // if the rmprimFile is not NULL, prepare to remove primitives
1743 } // remove primitives as desired by the user
1744
1745 FILE *fignore = fopen("IgnorePrims.info", "w");
1746 if (fignore == NULL) {
1747 printf(
1748 "error opening IgnorePrims.info file in write mode ... returning\n");
1749 return (-1);
1750 }
1751
1752 for (int prim = 1; prim <= OrgnlNbPrimitives; ++prim) {
1753 fprintf(fignore, "%d %d %d\n", OrgnlToEffPrim[prim][0],
1754 OrgnlToEffPrim[prim][1], OrgnlToEffPrim[prim][2]);
1755 }
1756
1757 fclose(fignore);
1758 } // Ignore unnecessary primitives from the final count
1759
1760 // Reduced-Order Modelling information
1761 printf("ROM: switch to primitive representation after %d repetitions.\n",
1762 PrimAfter);
1763
1764 // Store model data in native neBEM format
1765 char NativeInFile[256];
1766
1767 strcpy(NativeInFile, NativeOutDir);
1768 strcat(NativeInFile, "/neBEMNative.inp");
1769 FILE *fNativeInFile = fopen(NativeInFile, "w");
1770 fprintf(fNativeInFile, "#====>Input directory\n");
1771 fprintf(fNativeInFile, "%s\n", NativePrimDir);
1772 fprintf(fNativeInFile, "#====>No. of primitives:\n");
1773 fprintf(fNativeInFile, "%d\n", NbPrimitives);
1774 fprintf(fNativeInFile, "#====>No. of volumes:\n");
1775 fprintf(fNativeInFile, "%d\n", VolMax);
1776 for (int prim = 1; prim <= NbPrimitives; ++prim) {
1777 char NativePrimFile[256];
1778 char strPrimNb[11];
1779 sprintf(strPrimNb, "%d", prim);
1780 strcpy(NativePrimFile, "Primitive");
1781 strcat(NativePrimFile, strPrimNb);
1782 strcat(NativePrimFile, ".inp");
1783
1784 fprintf(fNativeInFile, "#Input filename for primitive:\n");
1785 fprintf(fNativeInFile, "%s\n", NativePrimFile);
1786
1787 strcpy(NativePrimFile, NativePrimDir);
1788 strcat(NativePrimFile, "/Primitive");
1789 strcat(NativePrimFile, strPrimNb);
1790 strcat(NativePrimFile, ".inp");
1791
1792 FILE *fNativePrim = fopen(NativePrimFile, "w");
1793
1794 fprintf(fNativePrim, "#Number of vertices:\n");
1795 fprintf(fNativePrim, "%d\n", NbVertices[prim]);
1796 for (int vertex = 0; vertex < NbVertices[prim]; ++vertex) {
1797 fprintf(fNativePrim, "#Vertex %d:\n", vertex);
1798 fprintf(fNativePrim, "%lg %lg %lg\n", XVertex[prim][vertex],
1799 YVertex[prim][vertex], ZVertex[prim][vertex]);
1800 } // for vertex
1801 fprintf(fNativePrim, "#Outward normal:\n");
1802 fprintf(fNativePrim, "%lg %lg %lg\n", XNorm[prim], YNorm[prim],
1803 ZNorm[prim]);
1804 fprintf(fNativePrim, "#Volume references:\n");
1805 fprintf(fNativePrim, "%d %d\n", VolRef1[prim], VolRef2[prim]);
1806 fprintf(fNativePrim, "#Maximum number of segments:\n");
1807 fprintf(fNativePrim, "%d %d\n", 0, 0); // use the trio target, min, max
1808 fprintf(fNativePrim, "#Information on X periodicity:\n");
1809 fprintf(fNativePrim, "%d %d %lg\n", PeriodicTypeX[prim], PeriodicInX[prim],
1810 XPeriod[prim]);
1811 fprintf(fNativePrim, "#Information on Y periodicity:\n");
1812 fprintf(fNativePrim, "%d %d %lg\n", PeriodicTypeY[prim], PeriodicInY[prim],
1813 YPeriod[prim]);
1814 fprintf(fNativePrim, "#Information on Z periodicity:\n");
1815 fprintf(fNativePrim, "%d %d %lg\n", PeriodicTypeZ[prim], PeriodicInZ[prim],
1816 ZPeriod[prim]);
1817
1818 fclose(fNativePrim);
1819 } // for prim
1820
1821 fprintf(fNativeInFile, "#====>No. of boundary conditions per element:\n");
1822 fprintf(fNativeInFile, "%d\n", 1); // CHECK!!!
1823 fprintf(fNativeInFile, "#====>Device output directory name:\n");
1824 fprintf(fNativeInFile, "NativeResults\n");
1825 fprintf(fNativeInFile, "#====>Map input file:\n");
1826 fprintf(fNativeInFile, "MapFile.inp\n");
1827 fclose(fNativeInFile);
1828
1829 for (int volume = 0; volume <= VolMax; ++volume) {
1830 // Note that materials from 1 to 10 are conductors and
1831 // from 11 to 20 are dielectrics
1832 int shape, material, boundarytype;
1833 double epsilon, potential, charge;
1834 neBEMVolumeDescription(volume, &shape, &material, &epsilon, &potential,
1835 &charge, &boundarytype);
1836
1837 char NativeVolFile[256];
1838 char strVolNb[11];
1839 sprintf(strVolNb, "%d", volume);
1840
1841 strcpy(NativeVolFile, NativePrimDir);
1842 strcat(NativeVolFile, "/Volume");
1843 strcat(NativeVolFile, strVolNb);
1844 strcat(NativeVolFile, ".inp");
1845
1846 FILE *fNativeVol = fopen(NativeVolFile, "w");
1847
1848 fprintf(fNativeVol, "#Shape of the volume:\n");
1849 fprintf(fNativeVol, "%d\n", shape);
1850 fprintf(fNativeVol, "#Material of the volume:\n");
1851 fprintf(fNativeVol, "%d\n", material);
1852 fprintf(fNativeVol, "#Relative permittivity:\n");
1853 fprintf(fNativeVol, "%lg\n", epsilon);
1854 fprintf(fNativeVol, "#Potential:\n");
1855 fprintf(fNativeVol, "%lg\n", potential);
1856 fprintf(fNativeVol, "#Charge:\n");
1857 fprintf(fNativeVol, "%lg\n", charge);
1858 fprintf(fNativeVol, "#Boundary type:\n");
1859 fprintf(fNativeVol, "%d\n", boundarytype);
1860
1861 fclose(fNativeVol);
1862 }
1863
1864 // create a dummy map file
1865 {
1866 char NativeMapFile[256];
1867
1868 strcpy(NativeMapFile, NativePrimDir);
1869 strcat(NativeMapFile, "/MapFile.inp");
1870
1871 FILE *fNativeMap = fopen(NativeMapFile, "w");
1872
1873 fprintf(fNativeMap, "#Number of maps\n");
1874 fprintf(fNativeMap, "%d\n", 0);
1875 fprintf(fNativeMap, "#Number of lines\n");
1876 fprintf(fNativeMap, "%d\n", 0);
1877 fprintf(fNativeMap, "#Number of points\n");
1878 fprintf(fNativeMap, "%d\n", 0);
1879
1880 fclose(fNativeMap);
1881 }
1882
1883 // Store primitive related data in a file for a new model, if opted for
1885 if (OptFormattedFile) {
1886 fstatus = WritePrimitives();
1887 if (fstatus) {
1888 neBEMMessage("neBEMReadGeometry - problem writing Primtives.\n");
1889 return -1;
1890 }
1891 } // formatted file
1892
1893 if (OptUnformattedFile) {
1895 "neBEMReadGeometry - unformatted write not inplemented yet.\n");
1896 return -1;
1897 } // unformatted file
1898 } // store primitives
1899
1900 printf("neBEMReadGeometry: Geometry read!\n");
1901
1902 neBEMState = 3; // info about primitives read in after initialization and Nbs
1903
1904 stopClock = clock();
1906 printf("to read geometry\n");
1907
1908 return (0); // Success!
1909} // neBEMReadGeometry ends
int WritePrimitives(void)
int ReadPrimitives(void)
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 neBEMVolumeDescription(int ivol, int *shape, int *material, double *epsilon, double *potential, double *charge, int *boundarytype)
INTFACEGLOBAL int neBEMGetNbPrimitives(void)
INTFACEGLOBAL int OptPrintVolumeDetails
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 neBEMGetMirror(int prim, int *ix, int *jx, double *sx, int *iy, int *jy, double *sy, int *iz, int *jz, double *sz)
INTFACEGLOBAL int OptPrintPrimaryDetails
INTFACEGLOBAL int neBEMGetPeriodicities(int prim, int *ix, int *jx, double *sx, int *iy, int *jy, double *sy, int *iz, int *jz, double *sz)
neBEMGLOBAL int * PeriodicInY
Definition: neBEM.h:79
neBEMGLOBAL int * volMaterial
Definition: neBEM.h:63
neBEMGLOBAL int * NbElmntsOnPrim
Definition: neBEM.h:94
neBEMGLOBAL int * MirrorTypeZ
Definition: neBEM.h:81
neBEMGLOBAL double * MirrorDistYFromOrigin
Definition: neBEM.h:82
neBEMGLOBAL int OptStorePrimitives
Definition: neBEM.h:45
neBEMGLOBAL int MaxNbVertices
Definition: neBEM.h:59
neBEMGLOBAL int * BndPlaneInYMax
Definition: neBEM.h:85
neBEMGLOBAL double * XBndPlaneInXMax
Definition: neBEM.h:87
neBEMGLOBAL double * YBndPlaneInYMin
Definition: neBEM.h:86
neBEMGLOBAL int ** OrgnlToEffPrim
Definition: neBEM.h:69
neBEMGLOBAL double * ZPeriod
Definition: neBEM.h:80
neBEMGLOBAL double * VBndPlaneInZMin
Definition: neBEM.h:88
neBEMGLOBAL int * BndPlaneInXMax
Definition: neBEM.h:85
neBEMGLOBAL double * YBndPlaneInYMax
Definition: neBEM.h:87
neBEMGLOBAL double * PrimOriginZ
Definition: neBEM.h:73
neBEMGLOBAL int * PeriodicTypeX
Definition: neBEM.h:78
neBEMGLOBAL double * VBndPlaneInYMax
Definition: neBEM.h:89
neBEMGLOBAL double * ZBndPlaneInZMax
Definition: neBEM.h:87
neBEMGLOBAL int * PeriodicTypeY
Definition: neBEM.h:78
neBEMGLOBAL double * Epsilon1
Definition: neBEM.h:75
neBEMGLOBAL double * PrimOriginY
Definition: neBEM.h:73
neBEMGLOBAL double * MirrorDistZFromOrigin
Definition: neBEM.h:83
neBEMGLOBAL DirnCosn3D * PrimDC
Definition: neBEM.h:74
neBEMGLOBAL int * BndPlaneInZMin
Definition: neBEM.h:84
neBEMGLOBAL double * VBndPlaneInXMin
Definition: neBEM.h:88
neBEMGLOBAL int NbFloatingConductors
Definition: neBEM.h:124
neBEMGLOBAL double * XPeriod
Definition: neBEM.h:80
neBEMGLOBAL double * volEpsilon
Definition: neBEM.h:64
neBEMGLOBAL int * volRef
Definition: neBEM.h:63
neBEMGLOBAL double * volPotential
Definition: neBEM.h:64
neBEMGLOBAL double * VBndPlaneInZMax
Definition: neBEM.h:89
neBEMGLOBAL double * AvAsgndChDen
Definition: neBEM.h:90
neBEMGLOBAL double * XBndPlaneInXMin
Definition: neBEM.h:86
neBEMGLOBAL int * BndPlaneInXMin
Definition: neBEM.h:84
neBEMGLOBAL int * MirrorTypeY
Definition: neBEM.h:81
neBEMGLOBAL double * AvChDen
Definition: neBEM.h:90
neBEMGLOBAL int * MirrorTypeX
Definition: neBEM.h:81
neBEMGLOBAL int * BndPlaneInZMax
Definition: neBEM.h:85
neBEMGLOBAL int * BndPlaneInYMin
Definition: neBEM.h:84
neBEMGLOBAL double * PrimLX
Definition: neBEM.h:72
neBEMGLOBAL double * PrimOriginX
Definition: neBEM.h:73
neBEMGLOBAL int * PeriodicInZ
Definition: neBEM.h:79
neBEMGLOBAL int * PeriodicInX
Definition: neBEM.h:79
neBEMGLOBAL double * ZBndPlaneInZMin
Definition: neBEM.h:86
neBEMGLOBAL double * VBndPlaneInYMin
Definition: neBEM.h:88
neBEMGLOBAL double * volCharge
Definition: neBEM.h:64
neBEMGLOBAL int OrgnlNbPrimitives
Definition: neBEM.h:58
neBEMGLOBAL int * volBoundaryType
Definition: neBEM.h:63
neBEMGLOBAL int * PeriodicTypeZ
Definition: neBEM.h:78
neBEMGLOBAL int VolMax
Definition: neBEM.h:77
neBEMGLOBAL double * YPeriod
Definition: neBEM.h:80
neBEMGLOBAL double * PrimLZ
Definition: neBEM.h:72
neBEMGLOBAL int * volShape
Definition: neBEM.h:63
neBEMGLOBAL double * MirrorDistXFromOrigin
Definition: neBEM.h:82
neBEMGLOBAL double * VBndPlaneInXMax
Definition: neBEM.h:89
neBEMGLOBAL double * Epsilon2
Definition: neBEM.h:75
double ** dmatrix(long nrl, long nrh, long ncl, long nch)
Definition: nrutil.c:97
int ** imatrix(long nrl, long nrh, long ncl, long nch)
Definition: nrutil.c:121

◆ neBEMSetDefaults()

INTFACEGLOBAL int neBEMSetDefaults ( void  )

◆ neBEMSolve()

INTFACEGLOBAL int neBEMSolve ( void  )

Definition at line 4117 of file neBEMInterface.c.

4117 {
4118 int dbgFn = 0;
4119
4120 clock_t startSolveClock = clock();
4121
4122 if (TimeStep < 1) {
4123 neBEMMessage("neBEMSolve - TimeStep cannot be less than one!;\n");
4124 neBEMMessage(" Please put TimeStep = 1 for static problems.\n");
4125 }
4126
4127 if ((neBEMState == 5) || (neBEMState == 8)) {
4128 if (neBEMState == 8) // neBEMState 8 must have inverted flag on
4129 { // so it must be a case looking for solution with a NewBC
4130 if (NewBC == 0) {
4131 neBEMMessage("neBEMSolve - NewBC zero for neBEMState = 8!");
4132 neBEMMessage(" - Nothing to be done ... returning.");
4133 return -1;
4134 }
4135 }
4136
4137 if (NewModel) { // effectively, NewMesh = NewBC = NewPP = 1;
4138 int fstatus = ComputeSolution();
4139 if (fstatus != 0) {
4140 neBEMMessage("neBEMSolve - NewModel");
4141 return -1;
4142 }
4143 } else { // NewModel == 0
4144 if (NewMesh) {
4145 // effectively, NewBC = NewPP = 1;
4146 int fstatus = ComputeSolution();
4147 if (fstatus != 0) {
4148 neBEMMessage("neBEMSolve - NewMesh");
4149 return -1;
4150 }
4151 } else { // NewModel == NewMesh == 0
4152 if (NewBC) { // effectively, NewPP = 1;
4153 int fstatus = ComputeSolution();
4154 if (fstatus != 0) {
4155 neBEMMessage("neBEMSolve - Failure computing new solution");
4156 return -1;
4157 }
4158 } else { // NewBC == 0
4159 if (NewPP) {
4160 int fstatus = ReadSolution();
4161 if (fstatus != 0) {
4162 neBEMMessage("neBEMSolve - Failure reading solution");
4163 return (-1);
4164 }
4165 } else { // NewPP == 0
4166 printf("neBEMSolve: Nothing to do ... returning ...\n");
4167 return (-1);
4168 } // NewPP == 0
4169 } // NewBC == 0
4170 } // NewModel == NewDiscretization == 0
4171 } // NewModel == 0
4172
4173 neBEMState = 9;
4174 } else {
4175 printf("neBEMSolve: neBEMSolve can be called only in state 5 / 8 ...\n");
4176 printf("returning ...\n");
4177 return (-1);
4178 }
4179
4180 if (FailureCntr) {
4181 printf(
4182 "neBEMSolve: Approximations were made while computing the influence "
4183 "coefficients.\n");
4184 printf(" Please check the \"%s/Isles.log\" file.\n", PPOutDir);
4185 }
4186
4187 clock_t stopSolveClock = clock();
4188 neBEMTimeElapsed(startSolveClock, stopSolveClock);
4189 printf("to complete solve\n");
4190
4191 // Prepare voxelized data that will be exported to Garfield++
4192 if (OptVoxel) {
4193 clock_t startVoxelClock = clock();
4194
4195 int fstatus = VoxelFPR();
4196 if (fstatus != 0) {
4197 neBEMMessage("neBEMSolve - Failure computing VoxelFPR");
4198 return -1;
4199 }
4200
4201 clock_t stopVoxelClock = clock();
4202 neBEMTimeElapsed(startVoxelClock, stopVoxelClock);
4203 printf("to compute VoxelFPR\n");
4204 }
4205
4206 // Prepare 3dMap data that will be exported to Garfield++
4207 if (OptMap) {
4208 clock_t startMapClock = clock();
4209
4210 int fstatus = MapFPR();
4211 if (fstatus != 0) {
4212 neBEMMessage("neBEMSolve - Failure computing MapFPR");
4213 return -1;
4214 }
4215
4216 clock_t stopMapClock = clock();
4217 neBEMTimeElapsed(startMapClock, stopMapClock);
4218 printf("to compute MapFPR\n");
4219 }
4220
4221 // allocate memory for potential and field components within FAST volume mesh
4222 // and compute / read FastVol data
4223 // Similar allocation, computation and reading may be necessary for the KnCh
4224 // effects.
4225 // The other approach could be to create fast volumes that always have both
4226 // the influences (elements + known charges) added together. This approach
4227 // seems more managable now and is being followed.
4228 // Thus, if we want to inspect the effects of elements and known charges
4229 // separately, we will have to generate one fast volume with OptKnCh = 0,
4230 // and another with OptKnCh = 1. Subtraction of these two fast volumes will
4231 // provide us with the effect of KnCh.
4232 if (OptFastVol) {
4233 int MaxXCells = BlkNbXCells[1];
4234 int MaxYCells = BlkNbYCells[1];
4235 int MaxZCells = BlkNbZCells[1];
4236 clock_t startFastClock = clock();
4237
4238 // find maximum number of Xcells etc in all the blocks
4239 // simplifies memory allocation using nrutils but hogs memory!
4240 for (int block = 1; block <= FastVol.NbBlocks; ++block) {
4241 if (block == 1) {
4242 MaxXCells = BlkNbXCells[1];
4243 MaxYCells = BlkNbYCells[1];
4244 MaxZCells = BlkNbZCells[1];
4245 } else {
4246 if (MaxXCells < BlkNbXCells[block]) MaxXCells = BlkNbXCells[block];
4247 if (MaxYCells < BlkNbYCells[block]) MaxYCells = BlkNbYCells[block];
4248 if (MaxZCells < BlkNbZCells[block]) MaxZCells = BlkNbZCells[block];
4249 }
4250 } // loop block for finding maxm cells among all the blocks
4251
4252 if (dbgFn) {
4253 printf("OptFastVol: %d\n", OptFastVol);
4254 printf("NbPtSkip: %d\n", NbPtSkip);
4255 printf("OptStaggerFastVol: %d\n", OptStaggerFastVol);
4256 printf("NbStgPtSkip: %d\n", NbStgPtSkip);
4257 printf("OptReadFastPF: %d\n", OptReadFastPF);
4258 printf("LX: %le\n", FastVol.LX);
4259 printf("LY: %le\n", FastVol.LY);
4260 printf("LZ: %le\n", FastVol.LZ);
4261 printf("CornerX: %le\n", FastVol.CrnrX);
4262 printf("CornerY: %le\n", FastVol.CrnrY);
4263 printf("CornerZ: %le\n", FastVol.CrnrZ);
4264 printf("YStagger: %le\n", FastVol.YStagger);
4265 printf("NbOfBlocks: %d\n", FastVol.NbBlocks);
4266 for (int block = 1; block <= FastVol.NbBlocks; ++block) {
4267 printf("NbOfXCells[%d]: %d\n", block, BlkNbXCells[block]);
4268 printf("NbOfYCells[%d]: %d\n", block, BlkNbYCells[block]);
4269 printf("NbOfZCells[%d]: %d\n", block, BlkNbZCells[block]);
4270 printf("LZ[%d]: %le\n", block, BlkLZ[block]);
4271 printf("CornerZ[%d]: %le\n", block, BlkCrnrZ[block]);
4272 }
4273 printf("NbOfOmitVols: %d\n", FastVol.NbOmitVols);
4274 if (FastVol.NbOmitVols) {
4275 for (int omit = 1; omit <= FastVol.NbOmitVols; ++omit) {
4276 printf("OmitVolLX[%d]: %le\n", omit, OmitVolLX[omit]);
4277 printf("OmitVolLY[%d]: %le\n", omit, OmitVolLY[omit]);
4278 printf("OmitVolLZ[%d]: %le\n", omit, OmitVolLZ[omit]);
4279 printf("OmitCrnrX[%d]: %le\n", omit, OmitVolCrnrX[omit]);
4280 printf("OmitCrnrY[%d]: %le\n", omit, OmitVolCrnrY[omit]);
4281 printf("OmitCrnrZ[%d]: %le\n", omit, OmitVolCrnrZ[omit]);
4282 }
4283 }
4284 printf("MaxXCells, MaxYCells, MaxZCells: %d, %d, %d\n", MaxXCells,
4285 MaxYCells, MaxZCells);
4286 } // dbgFn
4287
4288 /* Memory wastage!!! Improve as soon as possible. */
4289 FastPot = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1, MaxYCells + 1,
4290 1, MaxZCells + 1);
4291 FastFX = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1, MaxYCells + 1,
4292 1, MaxZCells + 1);
4293 FastFY = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1, MaxYCells + 1,
4294 1, MaxZCells + 1);
4295 FastFZ = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1, MaxYCells + 1,
4296 1, MaxZCells + 1);
4297
4298 if (OptStaggerFastVol) {
4299 /* Memory wastage!!! Improve as soon as possible. */
4300 FastStgPot = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1,
4301 MaxYCells + 1, 1, MaxZCells + 1);
4302 FastStgFX = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1,
4303 MaxYCells + 1, 1, MaxZCells + 1);
4304 FastStgFY = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1,
4305 MaxYCells + 1, 1, MaxZCells + 1);
4306 FastStgFZ = d4tensor(1, FastVol.NbBlocks, 1, MaxXCells + 1, 1,
4307 MaxYCells + 1, 1, MaxZCells + 1);
4308 } // if OptStaggerFastVol
4309
4310 if ((OptCreateFastPF) && (!OptReadFastPF)) // reading overrides creation
4311 {
4312 int fstatus = FastVolPF();
4313 if (fstatus != 0) {
4314 neBEMMessage("neBEMSolve - Failure computing FastVolPF");
4315 return -1;
4316 }
4317 } // if OptCreateFastPF
4318
4319 if (OptReadFastPF) {
4320 int nbXCells, nbYCells, nbZCells;
4321 int tmpblk;
4322 double xpt, ypt, zpt;
4323
4324 char FastVolPFFile[256];
4325 strcpy(FastVolPFFile, BCOutDir);
4326 strcat(FastVolPFFile, "/FastVolPF.out");
4327 FILE* fFastVolPF = fopen(FastVolPFFile, "r");
4328 if (fFastVolPF == NULL) {
4329 neBEMMessage("in neBEMSolve - FastVolPFFile");
4330 return -1;
4331 }
4332
4333 fscanf(fFastVolPF, "#block\tX\tY\tZ\tPot\tFX\tFY\tFZ\n");
4334
4335 for (int block = 1; block <= FastVol.NbBlocks; ++block) {
4336 nbXCells = BlkNbXCells[block];
4337 nbYCells = BlkNbYCells[block];
4338 nbZCells = BlkNbZCells[block];
4339
4340 for (int i = 1; i <= nbXCells + 1; ++i) {
4341 for (int j = 1; j <= nbYCells + 1; ++j) {
4342 for (int k = 1; k <= nbZCells + 1; ++k) {
4343 fscanf(fFastVolPF, "%d\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",
4344 &tmpblk, &xpt, &ypt, &zpt, &FastPot[block][i][j][k],
4345 &FastFX[block][i][j][k], &FastFY[block][i][j][k],
4346 &FastFZ[block][i][j][k]);
4347 } // loop k
4348 } // loop j
4349 } // loop i
4350 } // loop block
4351 fclose(fFastVolPF);
4352
4353 if (OptStaggerFastVol) {
4354 char FastStgVolPFFile[256];
4355 strcpy(FastStgVolPFFile, BCOutDir);
4356 strcat(FastStgVolPFFile, "/FastStgVolPF.out");
4357 FILE *fFastStgVolPF = fopen(FastStgVolPFFile, "r");
4358 if (fFastStgVolPF == NULL) {
4359 neBEMMessage("in neBEMSolve - FastStgVolPFFile");
4360 return -1;
4361 }
4362
4363 fscanf(fFastStgVolPF, "#block\tX\tY\tZ\tPot\tFX\tFY\tFZ\n");
4364
4365 for (int block = 1; block <= FastVol.NbBlocks; ++block) {
4366 nbXCells = BlkNbXCells[block];
4367 nbYCells = BlkNbYCells[block];
4368 nbZCells = BlkNbZCells[block];
4369
4370 for (int i = 1; i <= nbXCells + 1; ++i) {
4371 for (int j = 1; j <= nbYCells + 1; ++j) {
4372 for (int k = 1; k <= nbZCells + 1; ++k) {
4373 fscanf(fFastStgVolPF, "%d\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",
4374 &tmpblk, &xpt, &ypt, &zpt, &FastStgPot[block][i][j][k],
4375 &FastStgFX[block][i][j][k], &FastStgFY[block][i][j][k],
4376 &FastStgFZ[block][i][j][k]);
4377 } // loop k
4378 } // loop j
4379 } // loop i
4380 } // loop block
4381 fclose(fFastStgVolPF);
4382 } // if OptStaggerFastVol
4383 } // if OptReadFastPF
4384
4385 clock_t stopFastClock = clock();
4386 neBEMTimeElapsed(startFastClock, stopFastClock);
4387 printf("to compute / read FastVolPF\n");
4388 } // if OptFastVol
4389
4390 if (OptWtFldFastVol) // allocate memory for weighting field fast volume
4391 // variables
4392 {
4393 int MaxXCells = WtFldBlkNbXCells[1];
4394 int MaxYCells = WtFldBlkNbYCells[1];
4395 int MaxZCells = WtFldBlkNbZCells[1];
4396 clock_t startFastClock = clock();
4397
4398 // find maximum number of Xcells etc in all the blocks
4399 // simplifies memory allocation using nrutils but hogs memory!
4400 for (int block = 1; block <= WtFldFastVol.NbBlocks; ++block) {
4401 if (block == 1) {
4402 MaxXCells = WtFldBlkNbXCells[1];
4403 MaxYCells = WtFldBlkNbYCells[1];
4404 MaxZCells = WtFldBlkNbZCells[1];
4405 } else {
4406 if (MaxXCells < WtFldBlkNbXCells[block])
4407 MaxXCells = WtFldBlkNbXCells[block];
4408 if (MaxYCells < WtFldBlkNbYCells[block])
4409 MaxYCells = WtFldBlkNbYCells[block];
4410 if (MaxZCells < WtFldBlkNbZCells[block])
4411 MaxZCells = WtFldBlkNbZCells[block];
4412 }
4413 } // loop block for finding maxm cells among all the blocks
4414
4415 if (dbgFn) {
4416 printf("OptWtFldFastVol: %d\n", OptWtFldFastVol);
4417 printf("WtFldNbPtSkip: %d\n", WtFldNbPtSkip);
4418 printf("OptWtFldStaggerFastVol: %d\n", OptWtFldStaggerFastVol);
4419 printf("WtFldNbStgPtSkip: %d\n", WtFldNbStgPtSkip);
4420 printf("OptWtFldReadFastPF: %d\n", OptWtFldReadFastPF);
4421 printf("LX: %le\n", WtFldFastVol.LX);
4422 printf("LY: %le\n", WtFldFastVol.LY);
4423 printf("LZ: %le\n", WtFldFastVol.LZ);
4424 printf("CornerX: %le\n", WtFldFastVol.CrnrX);
4425 printf("CornerY: %le\n", WtFldFastVol.CrnrY);
4426 printf("CornerZ: %le\n", WtFldFastVol.CrnrZ);
4427 printf("YStagger: %le\n", WtFldFastVol.YStagger);
4428 printf("NbOfBlocks: %d\n", WtFldFastVol.NbBlocks);
4429 for (int block = 1; block <= WtFldFastVol.NbBlocks; ++block) {
4430 printf("NbOfXCells[%d]: %d\n", block, WtFldBlkNbXCells[block]);
4431 printf("NbOfYCells[%d]: %d\n", block, WtFldBlkNbYCells[block]);
4432 printf("NbOfZCells[%d]: %d\n", block, WtFldBlkNbZCells[block]);
4433 printf("LZ[%d]: %le\n", block, WtFldBlkLZ[block]);
4434 printf("CornerZ[%d]: %le\n", block, WtFldBlkCrnrZ[block]);
4435 }
4436 printf("NbOfOmitVols: %d\n", WtFldFastVol.NbOmitVols);
4438 for (int omit = 1; omit <= WtFldFastVol.NbOmitVols; ++omit) {
4439 printf("OmitVolLX[%d]: %le\n", omit, WtFldOmitVolLX[omit]);
4440 printf("OmitVolLY[%d]: %le\n", omit, WtFldOmitVolLY[omit]);
4441 printf("OmitVolLZ[%d]: %le\n", omit, WtFldOmitVolLZ[omit]);
4442 printf("OmitCrnrX[%d]: %le\n", omit, WtFldOmitVolCrnrX[omit]);
4443 printf("OmitCrnrY[%d]: %le\n", omit, WtFldOmitVolCrnrY[omit]);
4444 printf("OmitCrnrZ[%d]: %le\n", omit, WtFldOmitVolCrnrZ[omit]);
4445 }
4446 }
4447 printf("MaxXCells, MaxYCells, MaxZCells: %d, %d, %d\n", MaxXCells,
4448 MaxYCells, MaxZCells);
4449 } // dbgFn
4450
4451 /* Memory wastage!!! Improve as soon as possible. */
4452 WtFldFastPot = d4tensor(1, WtFldFastVol.NbBlocks, 1, MaxXCells + 1, 1,
4453 MaxYCells + 1, 1, MaxZCells + 1);
4454 WtFldFastFX = d4tensor(1, WtFldFastVol.NbBlocks, 1, MaxXCells + 1, 1,
4455 MaxYCells + 1, 1, MaxZCells + 1);
4456 WtFldFastFY = d4tensor(1, WtFldFastVol.NbBlocks, 1, MaxXCells + 1, 1,
4457 MaxYCells + 1, 1, MaxZCells + 1);
4458 WtFldFastFZ = d4tensor(1, WtFldFastVol.NbBlocks, 1, MaxXCells + 1, 1,
4459 MaxYCells + 1, 1, MaxZCells + 1);
4460
4462 /* Memory wastage!!! Improve as soon as possible. */
4463 WtFldFastStgPot = d4tensor(1, WtFldFastVol.NbBlocks, 1, MaxXCells + 1, 1,
4464 MaxYCells + 1, 1, MaxZCells + 1);
4465 WtFldFastStgFX = d4tensor(1, WtFldFastVol.NbBlocks, 1, MaxXCells + 1, 1,
4466 MaxYCells + 1, 1, MaxZCells + 1);
4467 WtFldFastStgFY = d4tensor(1, WtFldFastVol.NbBlocks, 1, MaxXCells + 1, 1,
4468 MaxYCells + 1, 1, MaxZCells + 1);
4469 WtFldFastStgFZ = d4tensor(1, WtFldFastVol.NbBlocks, 1, MaxXCells + 1, 1,
4470 MaxYCells + 1, 1, MaxZCells + 1);
4471 } // if OptWtFldStaggerFastVol
4472
4473 if ((OptWtFldCreateFastPF) && (!OptWtFldReadFastPF)) // reading overrides
4474 {
4475 // Computing weighting field fast volume has not been implemented
4477 "neBEMSolve - Failure computing WtFldFastVolPF: not implemented");
4478 return -1;
4479 } // if OptWtFldCreateFastPF
4480
4481 if (OptWtFldReadFastPF) // reading option overrides creation
4482 {
4483 int nbXCells, nbYCells, nbZCells;
4484 int tmpblk;
4485 double xpt, ypt, zpt;
4486
4487 char FastVolPFFile[256];
4488 strcpy(FastVolPFFile, BCOutDir);
4489 strcat(FastVolPFFile, "/WtFldFastVolPF.out");
4490 FILE *fFastVolPF = fopen(FastVolPFFile, "r");
4491 if (fFastVolPF == NULL) {
4492 neBEMMessage("in neBEMSolve - WtFldFastVolPFFile");
4493 return -1;
4494 }
4495
4496 fscanf(fFastVolPF, "#block\tX\tY\tZ\tPot\tFX\tFY\tFZ\n");
4497
4498 for (int block = 1; block <= WtFldFastVol.NbBlocks; ++block) {
4499 nbXCells = WtFldBlkNbXCells[block];
4500 nbYCells = WtFldBlkNbYCells[block];
4501 nbZCells = WtFldBlkNbZCells[block];
4502
4503 for (int i = 1; i <= nbXCells + 1; ++i) {
4504 for (int j = 1; j <= nbYCells + 1; ++j) {
4505 for (int k = 1; k <= nbZCells + 1; ++k) {
4506 fscanf(fFastVolPF, "%d\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",
4507 &tmpblk, &xpt, &ypt, &zpt, &WtFldFastPot[block][i][j][k],
4508 &WtFldFastFX[block][i][j][k], &WtFldFastFY[block][i][j][k],
4509 &WtFldFastFZ[block][i][j][k]);
4510 } // loop k
4511 } // loop j
4512 } // loop i
4513 } // loop block
4514 fclose(fFastVolPF);
4515
4517 char FastStgVolPFFile[256];
4518 FILE *fFastStgVolPF;
4519 strcpy(FastStgVolPFFile, BCOutDir);
4520 strcat(FastStgVolPFFile, "/WtFldFastStgVolPF.out");
4521 fFastStgVolPF = fopen(FastStgVolPFFile, "r");
4522
4523 if (fFastStgVolPF == NULL) {
4524 neBEMMessage("in neBEMSolve - WtFldFastStgVolPFFile");
4525 return -1;
4526 }
4527
4528 fscanf(fFastStgVolPF, "#block\tX\tY\tZ\tPot\tFX\tFY\tFZ\n");
4529
4530 for (int block = 1; block <= WtFldFastVol.NbBlocks; ++block) {
4531 nbXCells = WtFldBlkNbXCells[block];
4532 nbYCells = WtFldBlkNbYCells[block];
4533 nbZCells = WtFldBlkNbZCells[block];
4534
4535 for (int i = 1; i <= nbXCells + 1; ++i) {
4536 for (int j = 1; j <= nbYCells + 1; ++j) {
4537 for (int k = 1; k <= nbZCells + 1; ++k) {
4538 fscanf(fFastStgVolPF, "%d\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",
4539 &tmpblk, &xpt, &ypt, &zpt,
4540 &WtFldFastStgPot[block][i][j][k],
4541 &WtFldFastStgFX[block][i][j][k],
4542 &WtFldFastStgFY[block][i][j][k],
4543 &WtFldFastStgFZ[block][i][j][k]);
4544 } // loop k
4545 } // loop j
4546 } // loop i
4547 } // loop block
4548 fclose(fFastStgVolPF);
4549 } // if OptWtFldStaggerFastVol
4550 } // if OptWtFldReadFastPF
4551
4552 clock_t stopFastClock = clock();
4553 neBEMTimeElapsed(startFastClock, stopFastClock);
4554 printf("to compute / read FastVolPF\n");
4555 } // if OptWtFldFastVol
4556
4557 return (0);
4558} // neBEMSolve ends
int FastVolPF(void)
int VoxelFPR(void)
int MapFPR(void)
int ComputeSolution(void)
Definition: neBEM.c:36
int ReadSolution(void)
Definition: neBEM.c:3618
neBEMGLOBAL double **** FastFY
Definition: neBEM.h:444
neBEMGLOBAL double **** FastFZ
Definition: neBEM.h:444
neBEMGLOBAL double **** WtFldFastStgFZ
Definition: neBEM.h:515
neBEMGLOBAL double **** WtFldFastFZ
Definition: neBEM.h:513
neBEMGLOBAL double **** FastStgFZ
Definition: neBEM.h:446
neBEMGLOBAL double **** FastStgPot
Definition: neBEM.h:445
neBEMGLOBAL double **** WtFldFastFX
Definition: neBEM.h:513
neBEMGLOBAL double **** WtFldFastPot
Definition: neBEM.h:512
neBEMGLOBAL double **** FastStgFY
Definition: neBEM.h:446
neBEMGLOBAL double **** WtFldFastStgPot
Definition: neBEM.h:514
neBEMGLOBAL double **** FastFX
Definition: neBEM.h:444
neBEMGLOBAL double **** WtFldFastFY
Definition: neBEM.h:513
neBEMGLOBAL double **** WtFldFastStgFY
Definition: neBEM.h:515
neBEMGLOBAL double **** FastStgFX
Definition: neBEM.h:446
neBEMGLOBAL int TimeStep
Definition: neBEM.h:241
neBEMGLOBAL double **** WtFldFastStgFX
Definition: neBEM.h:515
neBEMGLOBAL double **** FastPot
Definition: neBEM.h:443
double **** d4tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh, long nwl, long nwh)
Definition: nrutil.c:260

◆ neBEMTimeElapsed()

INTFACEGLOBAL double neBEMTimeElapsed ( clock_t  t0,
clock_t  t1 
)

Definition at line 5306 of file neBEMInterface.c.

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

Referenced by ComputeSolution(), InvertMatrix(), neBEMBoundaryConditions(), neBEMDiscretize(), neBEMReadGeometry(), and neBEMSolve().

◆ neBEMVolumeCharge()

INTFACEGLOBAL double neBEMVolumeCharge ( int  volume)

Definition at line 4711 of file neBEMInterface.c.

4711 {
4712 // Initialise the sum
4713 double sumcharge = 0.0;
4714
4715 // Loop over all elements
4716 for (int elem = 1; elem <= NbElements; ++elem) {
4717 // Find the primitive number for the element
4718 int prim = (EleArr + elem - 1)->PrimitiveNb;
4719
4720 // Find out to which volume this belongs to
4721 int volref = VolRef1[prim];
4722
4723 // Skip the rest if the volume is not right
4724 if (volref + 1 != volume) {
4725 continue;
4726 }
4727
4728 // Calculate the periodic weight of the primitive
4729 int rptCnt = 0;
4730 if (PeriodicInX[prim] || PeriodicInY[prim] || PeriodicInZ[prim]) {
4731 for (int xrpt = -PeriodicInX[prim]; xrpt <= PeriodicInX[prim]; ++xrpt)
4732 for (int yrpt = -PeriodicInY[prim]; yrpt <= PeriodicInY[prim]; ++yrpt)
4733 for (int zrpt = -PeriodicInZ[prim]; zrpt <= PeriodicInZ[prim];
4734 ++zrpt) {
4735 if ((xrpt == 0) && (yrpt == 0) && (zrpt == 0))
4736 continue;
4737 else
4738 ++rptCnt;
4739 }
4740 } else {
4741 rptCnt = 1;
4742 }
4743
4744 // Add the charge
4745 // printf("Element: %d, volume: %d, charge: %g\n", elem, volref,
4746 // (EleArr+elem-1)->Solution * (EleArr+elem-1)->G.dA);
4747 sumcharge +=
4748 rptCnt * (EleArr + elem - 1)->Solution * (EleArr + elem - 1)->G.dA;
4749 } // loop over elements
4750
4751 // Return the result
4752 // printf("Charge on volume %d: %g C\n", volume, sumcharge);
4753 return sumcharge;
4754} // end of neBEMVolumeCharge
neBEMGLOBAL double * Solution
Definition: neBEM.h:237

◆ 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 
)

Referenced by MapFPR(), and VoxelFPR().

◆ neBEMVolumePrimitives()

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

◆ neBEMWeightingField()

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

Definition at line 4678 of file neBEMInterface.c.

4678 {
4679 double potential;
4680
4681 if (neBEMState < 9) {
4682 printf("neBEMWeightingField cannot be called before reaching state 9.\n");
4683 return (-1);
4684 }
4685
4686 if (OptFixedWtField) {
4687 // minimum computation, too restricted!
4688 potential = FixedWtPotential;
4689 field->X = FixedWtFieldX;
4690 field->Y = FixedWtFieldY;
4691 field->Z = FixedWtFieldZ;
4692 } else if (OptWtFldFastVol) {
4693 // bit more computation, lot more flexibility
4694 // Note: this is not the Creat or Read option
4695 int fstatus = WtFldFastPFAtPoint(point, &potential, field);
4696 if (fstatus != 0) {
4697 neBEMMessage("neBEMWeightingField - WtFldFastPFAtPoint");
4698 return DBL_MAX;
4699 }
4700 } else {
4701 int fstatus = WtPFAtPoint(point, &potential, field, IdWtField);
4702 if (fstatus != 0) {
4703 neBEMMessage("neBEMWeightingField - WtPFAtPoint");
4704 return DBL_MAX;
4705 }
4706 }
4707
4708 return potential;
4709} // neBEMWeightingField ends
int WtPFAtPoint(Point3D *globalP, double *Potential, Vector3D *globalF, int IdWtField)
int WtFldFastPFAtPoint(Point3D *globalP, double *Potential, Vector3D *globalF)

◆ ReadElements()

INTFACEGLOBAL int ReadElements ( void  )

Definition at line 5181 of file neBEMInterface.c.

5181 {
5182 char ElementFile[256];
5183
5184 strcpy(ElementFile, MeshOutDir);
5185 strcat(ElementFile, "/Elements/StoreElems.out");
5186
5187 FILE *fStrEle = fopen(ElementFile, "r");
5188 if (fStrEle == NULL) {
5189 neBEMMessage("ReadElements - Could not open file to read elements");
5190 return -1;
5191 }
5192
5193 fscanf(fStrEle, "%d %d\n", &NbSurfs, &NbWires);
5194
5195 for (int prim = 1; prim <= NbPrimitives; ++prim) {
5196 fscanf(fStrEle, "%d %d\n", &NbSurfSegX[prim], &NbSurfSegZ[prim]);
5197 fscanf(fStrEle, "%d\n", &NbWireSeg[prim]);
5198 }
5199
5200 fscanf(fStrEle, "%d\n", &NbElements);
5201
5202 if (neBEMState == 3) {
5203 if (EleArr) {
5204 Element *tmp = (Element *)realloc(EleArr, NbElements * sizeof(Element));
5205 if (tmp != NULL) {
5206 EleArr = tmp;
5207 EleCntr = 0;
5208 } else {
5209 free(EleArr);
5210 printf("neBEMDiscretize: Re-allocating EleArr failed.\n");
5211 return (1);
5212 }
5213 printf("neBEMDiscretize: Re-allocated EleArr.\n");
5214 } // if EleArr => re-allocation
5215 else {
5216 EleArr = (Element *)malloc(NbElements * sizeof(Element));
5217 if (EleArr == NULL) {
5218 neBEMMessage("neBEMDiscretize - EleArr malloc");
5219 return -1;
5220 }
5221 } // else EleArr => fresh allocation
5222 } else {
5223 neBEMMessage("neBEMDiscretize - EleArr malloc; neBEMState mismatch!");
5224 return -1;
5225 } // else neBEMState == 3
5226
5227 for (int ele = 1; ele <= NbElements; ++ele) {
5228 fscanf(fStrEle, "%hd %d %d %d %d\n", &(EleArr + ele - 1)->DeviceNb,
5229 &(EleArr + ele - 1)->ComponentNb, &(EleArr + ele - 1)->PrimitiveNb,
5230 &(EleArr + ele - 1)->InterfaceId, &(EleArr + ele - 1)->Id);
5231 fscanf(fStrEle, "%hd %le %le %le %le %le %le\n",
5232 &(EleArr + ele - 1)->G.Type, &(EleArr + ele - 1)->G.Origin.X,
5233 &(EleArr + ele - 1)->G.Origin.Y, &(EleArr + ele - 1)->G.Origin.Z,
5234 &(EleArr + ele - 1)->G.LX, &(EleArr + ele - 1)->G.LZ,
5235 &(EleArr + ele - 1)->G.dA);
5236 fscanf(fStrEle, "%le %le %le\n", &(EleArr + ele - 1)->G.DC.XUnit.X,
5237 &(EleArr + ele - 1)->G.DC.XUnit.Y,
5238 &(EleArr + ele - 1)->G.DC.XUnit.Z);
5239 fscanf(fStrEle, "%le %le %le\n", &(EleArr + ele - 1)->G.DC.YUnit.X,
5240 &(EleArr + ele - 1)->G.DC.YUnit.Y,
5241 &(EleArr + ele - 1)->G.DC.YUnit.Z);
5242 fscanf(fStrEle, "%le %le %le\n", &(EleArr + ele - 1)->G.DC.ZUnit.X,
5243 &(EleArr + ele - 1)->G.DC.ZUnit.Y,
5244 &(EleArr + ele - 1)->G.DC.ZUnit.Z);
5245 fscanf(fStrEle, "%hd %le\n", &(EleArr + ele - 1)->E.Type,
5246 &(EleArr + ele - 1)->E.Lambda);
5247 fscanf(fStrEle, "%hd %le %le %le %le\n", &(EleArr + ele - 1)->BC.NbOfBCs,
5248 &(EleArr + ele - 1)->BC.CollPt.X, &(EleArr + ele - 1)->BC.CollPt.Y,
5249 &(EleArr + ele - 1)->BC.CollPt.Z, &(EleArr + ele - 1)->BC.Value);
5250 fscanf(fStrEle, "%le %le\n", &(EleArr + ele - 1)->Solution,
5251 &(EleArr + ele - 1)->Assigned);
5252 }
5253
5254 fscanf(fStrEle, "%d %d %d %d\n", &NbPointsKnCh, &NbLinesKnCh, &NbAreasKnCh,
5255 &NbVolumesKnCh);
5256
5257 for (int pt = 1; pt <= NbPointsKnCh; ++pt) {
5258 fscanf(fStrEle, "%d %le\n", &(PointKnChArr + pt - 1)->Nb,
5259 &(PointKnChArr + pt - 1)->Assigned);
5260 fscanf(fStrEle, "%le %le %le\n", &(PointKnChArr + pt - 1)->P.X,
5261 &(PointKnChArr + pt - 1)->P.Y, &(PointKnChArr + pt - 1)->P.Z);
5262 }
5263
5264 for (int line = 1; line <= NbLinesKnCh; ++line) {
5265 fscanf(fStrEle, "%d %le %le\n", &(LineKnChArr + line - 1)->Nb,
5266 &(LineKnChArr + line - 1)->Radius,
5267 &(LineKnChArr + line - 1)->Assigned);
5268 fscanf(fStrEle, "%le %le %le\n", &(LineKnChArr + line - 1)->Start.X,
5269 &(LineKnChArr + line - 1)->Start.Y,
5270 &(LineKnChArr + line - 1)->Start.Z);
5271 fscanf(fStrEle, "%le %le %le\n", &(LineKnChArr + line - 1)->Stop.X,
5272 &(LineKnChArr + line - 1)->Stop.Y,
5273 &(LineKnChArr + line - 1)->Stop.Z);
5274 }
5275
5276 for (int area = 1; area <= NbAreasKnCh; ++area) {
5277 fscanf(fStrEle, "%d %d %le\n", &(AreaKnChArr + area - 1)->Nb,
5278 &(AreaKnChArr + area - 1)->NbVertices,
5279 &(AreaKnChArr + area - 1)->Assigned);
5280 for (int vert = 1; vert <= (AreaKnChArr + area - 1)->NbVertices; ++vert) {
5281 fscanf(fStrEle, "%le %le %le\n",
5282 &(AreaKnChArr + area - 1)->Vertex[vert].X,
5283 &(AreaKnChArr + area - 1)->Vertex[vert].Y,
5284 &(AreaKnChArr + area - 1)->Vertex[vert].Z);
5285 }
5286 }
5287
5288 for (int vol = 1; vol <= NbVolumesKnCh; ++vol) {
5289 fscanf(fStrEle, "%d %d %le\n", &(VolumeKnChArr + vol - 1)->Nb,
5290 &(VolumeKnChArr + vol - 1)->NbVertices,
5291 &(VolumeKnChArr + vol - 1)->Assigned);
5292 for (int vert = 1; vert <= (VolumeKnChArr + vol - 1)->NbVertices; ++vert) {
5293 fscanf(fStrEle, "%le %le %le\n",
5294 &(VolumeKnChArr + vol - 1)->Vertex[vert].X,
5295 &(VolumeKnChArr + vol - 1)->Vertex[vert].Y,
5296 &(VolumeKnChArr + vol - 1)->Vertex[vert].Z);
5297 }
5298 }
5299
5300 fclose(fStrEle);
5301
5302 return 0;
5303} // ReadElements ends

Referenced by neBEMDiscretize().

◆ ReadInitFile()

INTFACEGLOBAL int ReadInitFile ( char  filename[])

◆ ReadPrimitives()

INTFACEGLOBAL int ReadPrimitives ( void  )

Definition at line 5076 of file neBEMInterface.c.

5076 {
5077 int dbgFn = 0;
5078
5079 char PrimitiveFile[256];
5080
5081 strcpy(PrimitiveFile, ModelOutDir);
5082 strcat(PrimitiveFile, "/Primitives/StorePrims.out");
5083
5084 FILE *fStrPrm;
5085 fStrPrm = fopen(PrimitiveFile, "r");
5086 if (fStrPrm == NULL) {
5087 neBEMMessage("ReadPrimitives - Could not open file to read primitives");
5088 return -1;
5089 }
5090
5091 fscanf(fStrPrm, "%d %d\n", &NbVolumes, &VolMax);
5092 fscanf(fStrPrm, "%d\n", &NbPrimitives);
5093 fscanf(fStrPrm, "%d\n", &MaxNbVertices);
5094
5095 // assign neBEMState and allocate memory
5096 neBEMState = 2;
5105 Radius = dvector(1, NbPrimitives); // can lead to a little memory misuse
5110 NbWireSeg = ivector(1, NbPrimitives); // little memory misuse
5127
5128 for (int prim = 1; prim <= NbPrimitives; ++prim) {
5129 fscanf(fStrPrm, "%d\n", &PrimType[prim]);
5130 fscanf(fStrPrm, "%d\n", &InterfaceType[prim]);
5131 fscanf(fStrPrm, "%d\n", &NbVertices[prim]);
5132
5133 for (int vert = 0; vert < NbVertices[prim]; ++vert) {
5134 fscanf(fStrPrm, "%le %le %le\n", &XVertex[prim][vert],
5135 &YVertex[prim][vert], &ZVertex[prim][vert]);
5136 } // vert loop
5137
5138 fscanf(fStrPrm, "%le %le %le\n", &XNorm[prim], &YNorm[prim], &ZNorm[prim]);
5139 fscanf(fStrPrm, "%le\n", &Radius[prim]);
5140
5141 fscanf(fStrPrm, "%le %le %le %le %le\n", &Epsilon1[prim], &Epsilon2[prim],
5142 &Lambda[prim], &ApplPot[prim], &ApplCh[prim]);
5143
5144 fscanf(fStrPrm, "%d %d\n", &VolRef1[prim], &VolRef2[prim]);
5145
5146 fscanf(fStrPrm, "%d %d %d\n", &PeriodicTypeX[prim], &PeriodicTypeY[prim],
5147 &PeriodicTypeZ[prim]);
5148 fscanf(fStrPrm, "%d %d %d\n", &PeriodicInX[prim], &PeriodicInY[prim],
5149 &PeriodicInZ[prim]);
5150 fscanf(fStrPrm, "%le %le %le\n", &XPeriod[prim], &YPeriod[prim],
5151 &ZPeriod[prim]);
5152 fscanf(fStrPrm, "%le %le %le\n", &MirrorDistXFromOrigin[prim],
5154 } // prim loop
5155
5156 volRef = ivector(0, VolMax);
5157 volShape = ivector(0, VolMax);
5161 volCharge = dvector(0, VolMax);
5163 for (int volref = 0; volref <= VolMax; ++volref) {
5164 neBEMVolumeDescription(volref, &volShape[volref], &volMaterial[volref],
5165 &volEpsilon[volref], &volPotential[volref],
5166 &volCharge[volref], &volBoundaryType[volref]);
5167 if (dbgFn) {
5168 printf("volref: %d\n", volref);
5169 printf("shape: %d, material: %d\n", volShape[volref],
5170 volMaterial[volref]);
5171 printf("eps: %lg, pot: %lg\n", volEpsilon[volref], volPotential[volref]);
5172 printf("q: %lg, type: %d\n", volCharge[volref], volBoundaryType[volref]);
5173 }
5174 } // volume loop
5175
5176 fclose(fStrPrm);
5177
5178 return 0;
5179} // ReadPrimitives ends
neBEMGLOBAL int NbVolumes
Definition: neBEM.h:56

Referenced by neBEMReadGeometry().

◆ WriteElements()

INTFACEGLOBAL int WriteElements ( void  )

Definition at line 4982 of file neBEMInterface.c.

4982 {
4983 char ElementFile[256];
4984
4985 strcpy(ElementFile, MeshOutDir);
4986 strcat(ElementFile, "/Elements/StoreElems.out");
4987
4988 FILE *fStrEle = fopen(ElementFile, "w");
4989 if (fStrEle == NULL) {
4990 neBEMMessage("WriteElements - Could not create file to store elements");
4991 return -1;
4992 }
4993
4994 fprintf(fStrEle, "%d %d\n", NbSurfs, NbWires);
4995
4996 for (int prim = 1; prim <= NbPrimitives; ++prim) {
4997 fprintf(fStrEle, "%d %d\n", NbSurfSegX[prim], NbSurfSegZ[prim]);
4998 fprintf(fStrEle, "%d\n", NbWireSeg[prim]);
4999 }
5000
5001 fprintf(fStrEle, "%d\n", NbElements);
5002
5003 for (int ele = 1; ele <= NbElements; ++ele) {
5004 fprintf(fStrEle, "%d %d %d %d %d\n", (EleArr + ele - 1)->DeviceNb,
5005 (EleArr + ele - 1)->ComponentNb, (EleArr + ele - 1)->PrimitiveNb,
5006 (EleArr + ele - 1)->InterfaceId, (EleArr + ele - 1)->Id);
5007 fprintf(fStrEle, "%d %le %le %le %le %le %le\n", (EleArr + ele - 1)->G.Type,
5008 (EleArr + ele - 1)->G.Origin.X, (EleArr + ele - 1)->G.Origin.Y,
5009 (EleArr + ele - 1)->G.Origin.Z, (EleArr + ele - 1)->G.LX,
5010 (EleArr + ele - 1)->G.LZ, (EleArr + ele - 1)->G.dA);
5011 fprintf(fStrEle, "%le %le %le\n", (EleArr + ele - 1)->G.DC.XUnit.X,
5012 (EleArr + ele - 1)->G.DC.XUnit.Y, (EleArr + ele - 1)->G.DC.XUnit.Z);
5013 fprintf(fStrEle, "%le %le %le\n", (EleArr + ele - 1)->G.DC.YUnit.X,
5014 (EleArr + ele - 1)->G.DC.YUnit.Y, (EleArr + ele - 1)->G.DC.YUnit.Z);
5015 fprintf(fStrEle, "%le %le %le\n", (EleArr + ele - 1)->G.DC.ZUnit.X,
5016 (EleArr + ele - 1)->G.DC.ZUnit.Y, (EleArr + ele - 1)->G.DC.ZUnit.Z);
5017 fprintf(fStrEle, "%d %le\n", (EleArr + ele - 1)->E.Type,
5018 (EleArr + ele - 1)->E.Lambda);
5019 fprintf(fStrEle, "%d %le %le %le %le\n", (EleArr + ele - 1)->BC.NbOfBCs,
5020 (EleArr + ele - 1)->BC.CollPt.X, (EleArr + ele - 1)->BC.CollPt.Y,
5021 (EleArr + ele - 1)->BC.CollPt.Z, (EleArr + ele - 1)->BC.Value);
5022 fprintf(fStrEle, "%le %le\n", (EleArr + ele - 1)->Solution,
5023 (EleArr + ele - 1)->Assigned);
5024 }
5025
5026 fprintf(fStrEle, "%d %d %d %d\n", NbPointsKnCh, NbLinesKnCh, NbAreasKnCh,
5028
5029 for (int pt = 1; pt <= NbPointsKnCh; ++pt) {
5030 fprintf(fStrEle, "%d %le\n", (PointKnChArr + pt - 1)->Nb,
5031 (PointKnChArr + pt - 1)->Assigned);
5032 fprintf(fStrEle, "%le %le %le\n", (PointKnChArr + pt - 1)->P.X,
5033 (PointKnChArr + pt - 1)->P.Y, (PointKnChArr + pt - 1)->P.Z);
5034 }
5035
5036 for (int line = 1; line <= NbLinesKnCh; ++line) {
5037 fprintf(fStrEle, "%d %le %le\n", (LineKnChArr + line - 1)->Nb,
5038 (LineKnChArr + line - 1)->Radius,
5039 (LineKnChArr + line - 1)->Assigned);
5040 fprintf(fStrEle, "%le %le %le\n", (LineKnChArr + line - 1)->Start.X,
5041 (LineKnChArr + line - 1)->Start.Y,
5042 (LineKnChArr + line - 1)->Start.Z);
5043 fprintf(fStrEle, "%le %le %le\n", (LineKnChArr + line - 1)->Stop.X,
5044 (LineKnChArr + line - 1)->Stop.Y, (LineKnChArr + line - 1)->Stop.Z);
5045 }
5046
5047 for (int area = 1; area <= NbAreasKnCh; ++area) {
5048 fprintf(fStrEle, "%d %d %le\n", (AreaKnChArr + area - 1)->Nb,
5049 (AreaKnChArr + area - 1)->NbVertices,
5050 (AreaKnChArr + area - 1)->Assigned);
5051 for (int vert = 1; vert <= (AreaKnChArr + area - 1)->NbVertices; ++vert) {
5052 fprintf(fStrEle, "%le %le %le\n",
5053 (AreaKnChArr + area - 1)->Vertex[vert].X,
5054 (AreaKnChArr + area - 1)->Vertex[vert].Y,
5055 (AreaKnChArr + area - 1)->Vertex[vert].Z);
5056 }
5057 }
5058
5059 for (int vol = 1; vol <= NbVolumesKnCh; ++vol) {
5060 fprintf(fStrEle, "%d %d %le\n", (VolumeKnChArr + vol - 1)->Nb,
5061 (VolumeKnChArr + vol - 1)->NbVertices,
5062 (VolumeKnChArr + vol - 1)->Assigned);
5063 for (int vert = 1; vert <= (VolumeKnChArr + vol - 1)->NbVertices; ++vert) {
5064 fprintf(fStrEle, "%le %le %le\n",
5065 (VolumeKnChArr + vol - 1)->Vertex[vert].X,
5066 (VolumeKnChArr + vol - 1)->Vertex[vert].Y,
5067 (VolumeKnChArr + vol - 1)->Vertex[vert].Z);
5068 }
5069 }
5070
5071 fclose(fStrEle);
5072
5073 return 0;
5074} // WriteElements ends

Referenced by neBEMDiscretize().

◆ WritePrimitives()

INTFACEGLOBAL int WritePrimitives ( void  )

Definition at line 4933 of file neBEMInterface.c.

4933 {
4934 char PrimitiveFile[256];
4935
4936 strcpy(PrimitiveFile, ModelOutDir);
4937 strcat(PrimitiveFile, "/Primitives/StorePrims.out");
4938
4939 FILE *fStrPrm = fopen(PrimitiveFile, "w");
4940 if (fStrPrm == NULL) {
4941 neBEMMessage("WritePrimitives - Could not create file to store primitives");
4942 return -1;
4943 }
4944
4945 fprintf(fStrPrm, "%d %d\n", NbVolumes, VolMax);
4946 fprintf(fStrPrm, "%d\n", NbPrimitives);
4947 fprintf(fStrPrm, "%d\n", MaxNbVertices);
4948
4949 for (int prim = 1; prim <= NbPrimitives; ++prim) {
4950 fprintf(fStrPrm, "%d\n", PrimType[prim]);
4951 fprintf(fStrPrm, "%d\n", InterfaceType[prim]);
4952 fprintf(fStrPrm, "%d\n", NbVertices[prim]);
4953
4954 for (int vert = 0; vert < NbVertices[prim]; ++vert) {
4955 fprintf(fStrPrm, "%le %le %le\n", XVertex[prim][vert],
4956 YVertex[prim][vert], ZVertex[prim][vert]);
4957 } // vert loop
4958
4959 fprintf(fStrPrm, "%le %le %le\n", XNorm[prim], YNorm[prim], ZNorm[prim]);
4960 fprintf(fStrPrm, "%le\n", Radius[prim]);
4961
4962 fprintf(fStrPrm, "%le %le %le %le %le\n", Epsilon1[prim], Epsilon2[prim],
4963 Lambda[prim], ApplPot[prim], ApplCh[prim]);
4964
4965 fprintf(fStrPrm, "%d %d\n", VolRef1[prim], VolRef2[prim]);
4966
4967 fprintf(fStrPrm, "%d %d %d\n", PeriodicTypeX[prim], PeriodicTypeY[prim],
4968 PeriodicTypeZ[prim]);
4969 fprintf(fStrPrm, "%d %d %d\n", PeriodicInX[prim], PeriodicInY[prim],
4970 PeriodicInZ[prim]);
4971 fprintf(fStrPrm, "%le %le %le\n", XPeriod[prim], YPeriod[prim],
4972 ZPeriod[prim]);
4973 fprintf(fStrPrm, "%le %le %le\n", MirrorDistXFromOrigin[prim],
4975 } // prim loop
4976
4977 fclose(fStrPrm);
4978
4979 return 0;
4980} // 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 neBEMDiscretize().

◆ fgnuMesh

INTFACEGLOBAL FILE * fgnuMesh

Definition at line 55 of file neBEMInterface.h.

Referenced by neBEMDiscretize().

◆ fgnuPrim

INTFACEGLOBAL FILE* fgnuPrim

◆ NbThreads

INTFACEGLOBAL int NbThreads

Definition at line 49 of file neBEMInterface.h.

Referenced by neBEMInitialize().

◆ neBEMFieldCallCntr

INTFACEGLOBAL int neBEMFieldCallCntr

Definition at line 117 of file neBEMInterface.h.

◆ neBEMState

◆ OptDeviceFile

INTFACEGLOBAL int OptDeviceFile

Definition at line 44 of file neBEMInterface.h.

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

◆ OptElementFiles

◆ OptGnuplot

◆ OptGnuplotElements

◆ OptGnuplotPrimitives

◆ OptPrimitiveFiles

◆ OptPrintPrimaryDetails

INTFACEGLOBAL int OptPrintPrimaryDetails

◆ OptPrintVertexAndNormal

◆ OptPrintVolumeDetails

INTFACEGLOBAL int OptPrintVolumeDetails

◆ OptReuseDir

INTFACEGLOBAL int OptReuseDir

Definition at line 150 of file neBEMInterface.h.

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

◆ startClock

INTFACEGLOBAL clock_t startClock

◆ stopClock

INTFACEGLOBAL clock_t stopClock