int main(void)
{
neBEMInitialize();
neBEMReadGeometry();
neBEMDiscretize(elementNbs);
neBEMBoundaryConditions();
neBEMSolve();
PostProcess();
neBEMEnd();
}
Let us see one of the example codes (named ExampleDevice) to check the details. A brief discussion follows the source code which itself is self-explanatory.
/* Source code begins */
// Analyze an Example Device
#include <stdio.h>
#include <assert.h>
#include <Interface.h>
#include <Vector.h>
#include <NR.h>
#include <neBEM.h>
int main(void)
{
int PostProcess(void); // a typical post-processing function
neBEMState = 0;
int fstatus = neBEMInitialize();
assert(fstatus == 0);
fstatus = neBEMReadGeometry();
assert(fstatus == 0);
int **elementNbs;
elementNbs = imatrix(1, NbPrimitives, 1, 2);
for(register int prim = 1; prim <= NbPrimitives; ++prim)
{
elementNbs[prim][1] = tmpNbXSegs[prim];
elementNbs[prim][2] = tmpNbZSegs[prim];
}
fstatus = neBEMDiscretize(elementNbs);
assert(fstatus == 0);
fstatus = neBEMBoundaryConditions();
assert (fstatus == 0);
NewModel = 1; NewMesh = 1; NewBoundaryCondn = 1; NewPostProcess = 1;
ModelCntr = 1; MeshCntr = 1; BoundaryCondnCntr = 1; PostProcessCntr = 1;
TimeStep = 1; EndOfTime = 1;
fstatus = neBEMSolve(); // fresh calculation
assert(fstatus == 0);
fstatus = PostProcess();
assert(fstatus == 0);
fstatus = neBEMEnd();
assert(fstatus == 0);
return 0;
} // main ends
/* Source code ends */
At the beginning of neBEM computation, the user needs to assign the state variable neBEMState.
neBEMState = 0; int fstatus = neBEMInitialize(); assert(fstatus == 0);Within the neBEMInitialize function, default values of global variables used by neBEM are set using the
neBEMSetDefaultsfunction. If the OptDeviceFile variable is set to a non-zero value via neBEMSetDefaults, it is decided that the device details are to be read in from a file of name as assigned to a string named DeviceInputFile. In such an event neBEMInitialize invokes
neBEMGetInputFromFilefunction after setting up the defaults. For an example where the DeviceInputFile represents a real input file, check the examples related to ExampleDevice, microMegas and Iarocci tube. On the other hand, if OptDeviceFile is set to zero through neBEMSetDefaults, it is assumed that the user will supply all the necessary details related to the device model using hard-coded numbers, evaluation of suitable expressions within the interface source code, or through the execution of other function calls presumably residing in external packages. The GarfieldInterface, for example, follows the latter and invokes Garfield functions such as bemnpr, bempri etc to specify necessary details from the functions discussed below. This could be true for any other code being interfaced to the neBEM toolkit. For very simple examples that do not use external files for getting in the device details, check the examples related to Plate, DielectricInterface (DI, DItri etc).
The neBEMReadGeometry function uses neBEMGetNbPrimitives, neBEMGetPrimitive,
neBEMVolumeDescription and neBEMVolumePoint, the last one used sparingly, to
define the device details.
Function neBEMDisctrize sets up a mesh on the primitives used to define the device. It is likely that a routine will be supplied either by neBEM or the user to analyze the primitives and to decide the size and shape of elements necessary on each primitive. At present, the corase-ness or fine-ness is dependent on these user-inputs being supplied to neBEMDiscretize as elementNbs array. For a wire, the first element of this array needs to be 1 or more to be considered for meshing. For a surface, both the elements need to be 1 or more. If the above conditions for the elements of elementNbs array fail, the control on meshing goes to the first three variables specified in the init files, namely MinNbElementsOnLength, MaxNbElementsOnLength and ElementLengthRqstd, or the NEBEM procedure call in Garfield. We understand that while this is a convenient approach, the resulting mesh is unlikely to have good qualities. While it is true that neBEM can produce reasonably good results even on very coarse and uneven mesh, efforts will be made to make the toolkit intelligent enough to take correct decisions on discretization matters.
Since a boundary condition is likely to be associated with each element, the corresponding routine, neBEMBoundaryConditions can be called only after the primitives have been discretized.
Finally, the solution is requested after setting up several important flags that define the problem. These flags follow from the assumption that any unique device can have several models, each model can be discretized to different extent, each discretized device can be subjected to different boundary conditions leading to different singularity distributions, and, finally, each solution may be used to estimate various properties of a device configuration of interest. For example, a micromegas mesh in a micromegas TPC under study, can be modeled using surfaces or by using wires resulting into completely different sets of primitves. These primitives can have different meshing in order to attain different levels of accuracy. Till this far, the problem is determined almost entirely by the geometry of the problem and the material properties such as the dielectric permittivity of the components while analyzing an electrostatic problem. The discretized device can be subject to different boundary conditions each of which will result into distinct solutions, i.e., singularity distribution. These singularity distributions can be used to estimate different properties of the configuration under study by carrying out different post-processes. These flags have a hieracrchy, the first mentioned having the highest priority:
| Flag | Explanation |
|---|---|
| NewModel | 1 implies a fresh calculation. |
| NewMesh | 1 implies a new mesh for a device. |
| NewBoundaryCondn | 1 implies new RHS for the same LHS; |
| skips matrix inversion. | |
| NewPostProcess | 1 implies the use of the same solution; |
| skips matrix inversion, as well as | |
| the step for computing the solution. |
We maintain four counters as well:
| Counter | Explanation |
|---|---|
| ModelCntr | keeps track of the model for a given device |
| MeshCntr | keeps track of the mesh for a given model |
| BCCntr | keeps track of the association of the boundary |
| condition and its solution. This has to | |
| maintained by the user manually and supplied, | |
| for example, while carrying out a post-processing | |
| for a solution that was computed before. | |
| PPCntr | numbers different post-processes for a given |
| solution resulting from a given set of preceding | |
| conditions. |
If NewModel is set to 1, all the other flags are automatically ignored and a fresh calculation is initiated. The counters are, however, used only to create a new directory structure or use an existing one.
Please note that for solving a static problem, TimeStep and EndOfTime should both be put equal to one.
Post-processing can be carried out with separate function(s) written within this code, or function(s) residing in a separate source code, or a mix. Our experience indicate that the post-processing functions are quite complex entities themselves and are better off if kept independent of the driver routine. Among many things to be carried out in the post-processing phase, one could be the evaluation of weighting field. This is an important issue for nuclear detectors and can be carried out in two simple steps, as follows:
int IdWtField = neBEMPrepareWeightingField(int nprim, int primlist[]); int status = neBEMWeightingField(Point3D point, Vector3D *field, int IdWtField);
Here, IdWtField identifies the weighting field identification tag for which we have nprim number of primitives identified by the list primlist raised to 1, while all the other primitives are maintained grounded. At any point, the weighting field for the configuration identified by the tag IdWtField is obtained through the varaible field.
Finally, neBEMEnd ends neBEM gracefully, writing and closing necessary log files.