Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Macros | Functions
ReTriM.c File Reference
#include <stdio.h>
#include <time.h>
#include <unistd.h>
#include "neBEMInterface.h"
#include "Isles.h"
#include "NR.h"
#include "Vector.h"
#include "neBEM.h"

Go to the source code of this file.

Macros

#define MyPI   3.14159265358979323846
 

Functions

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)
 
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 NbSegs)
 
int AnalyzePrimitive (int prim, int *NbSegCoord1, int *NbSegCoord2)
 
int AnalyzeWire (int prim, int *NbSeg)
 
int AnalyzeSurface (int prim, int *NbSegCoord1, int *NbSegCoord2)
 
int DiscretizeWire (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 NbSegs)
 
int DiscretizeTriangle (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)
 
int DiscretizeRectangle (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)
 
int DiscretizePolygon (int, int, double[], double[], double[], double, double, double, int, int, int, double, double, double, int NbSegX, int NbSegZ)
 
int BoundaryConditions (void)
 

Macro Definition Documentation

◆ MyPI

#define MyPI   3.14159265358979323846

Definition at line 15 of file ReTriM.c.

Function Documentation

◆ AnalyzePrimitive()

int AnalyzePrimitive ( int  prim,
int *  NbSegCoord1,
int *  NbSegCoord2 
)

Definition at line 121 of file ReTriM.c.

121 {
122 int fstatus;
123
124 switch (NbVertices[prim]) {
125 case 2: // wire
126 fstatus = AnalyzeWire(prim, NbSegCoord1);
127 *NbSegCoord2 = 0;
128 // assert(fstatus == 0);
129 if (fstatus != 0) {
130 neBEMMessage("AnalyzePrimitive - AnalyzeWire");
131 return -1;
132 }
133 return (2);
134 break;
135 case 3: // triangle
136 fstatus = AnalyzeSurface(prim, NbSegCoord1, NbSegCoord2);
137 // assert(fstatus == 0);
138 if (fstatus != 0) {
139 neBEMMessage("AnalyzePrimitive - AnalyzeSurface");
140 return -1;
141 }
142 return (3);
143 break;
144 case 4: // rectangle
145 fstatus = AnalyzeSurface(prim, NbSegCoord1, NbSegCoord2);
146 // assert(fstatus == 0);
147 if (fstatus != 0) {
148 neBEMMessage("AnalyzePrimitive - AnalyzeSurface");
149 return -1;
150 }
151 return (4);
152 break;
153 default: // no shape!
154 return (0);
155 }
156} // end of AnalyzePrimitive
AnalyzeSurface
int AnalyzeSurface(int prim, int *NbSegCoord1, int *NbSegCoord2)
Definition: ReTriM.c:246
AnalyzeWire
int AnalyzeWire(int prim, int *NbSeg)
Definition: ReTriM.c:158
neBEMMessage
int neBEMMessage(const char *message)
Definition: neBEMInterface.c:4927
NbVertices
neBEMGLOBAL int * NbVertices
Definition: neBEM.h:70

Referenced by neBEMDiscretize().

◆ AnalyzeSurface()

int AnalyzeSurface ( int  prim,
int *  NbSegCoord1,
int *  NbSegCoord2 
)

Definition at line 246 of file ReTriM.c.

246 {
247 int nb1 = *NbSegCoord1, nb2 = *NbSegCoord2;
248
249 if ((nb1 < 1) || (nb2 < 1)) // absurd! use the trio: target, min, max
250 {
251 // Triangle primitives have their right angle on vertex 1
252 double l1 = (XVertex[prim][0] - XVertex[prim][1]) *
253 (XVertex[prim][0] - XVertex[prim][1]) +
254 (YVertex[prim][0] - YVertex[prim][1]) *
255 (YVertex[prim][0] - YVertex[prim][1]) +
256 (ZVertex[prim][0] - ZVertex[prim][1]) *
257 (ZVertex[prim][0] - ZVertex[prim][1]);
258 l1 = sqrt(l1);
259 double l2 = (XVertex[prim][2] - XVertex[prim][1]) *
260 (XVertex[prim][2] - XVertex[prim][1]) +
261 (YVertex[prim][2] - YVertex[prim][1]) *
262 (YVertex[prim][2] - YVertex[prim][1]) +
263 (ZVertex[prim][2] - ZVertex[prim][1]) *
264 (ZVertex[prim][2] - ZVertex[prim][1]);
265 l2 = sqrt(l2);
266
267 // We can use the lengths independently and forget about area
268 // for the time being
269
270 // double area = l1 * l2;
271
272 nb1 = (int)(l1 / ElementLengthRqstd);
273 if ((nb1 > MinNbElementsOnLength) &&
274 (nb1 < MaxNbElementsOnLength)) {
275 // nothing to be done
276 } else if (l1 < MINDIST) {
277 fprintf(fMeshLog, "Length1 too small on primitive %d!\n", prim);
278 nb1 = 1;
279 } else if (nb1 < MinNbElementsOnLength) {
280 // Need to have at least MinNbElementsOnLength elements per wire primitive
281 nb1 = MinNbElementsOnLength;
282 double ellength = l1 / (double)nb1;
283 // which may not be possible if the length is small
284 if (ellength < MINDIST) {
285 nb1 = (int)(l1 / MINDIST);
286 // However, it is necessary to have at least one element!
287 if (nb1 < 1) {
288 fprintf(fMeshLog, "Length1 very small on primitive %d!\n", prim);
289 nb1 = 1;
290 }
291 }
292 } // else if nb1 < MinNbElementsOnLength
293
294 if (nb1 > MaxNbElementsOnLength) {
295 fprintf(fMeshLog, "Too many elements on Length1 for primitive %d!\n",
296 prim);
297 fprintf(fMeshLog, "Number of elements reduced to maximum allowed %d\n",
298 MaxNbElementsOnLength);
299 nb1 = MaxNbElementsOnLength;
300 }
301
302 nb2 = (int)(l2 / ElementLengthRqstd);
303 if ((nb2 > MinNbElementsOnLength) &&
304 (nb2 < MaxNbElementsOnLength)) {
305 // nothing to be done
306 } else if (l2 < MINDIST) {
307 fprintf(fMeshLog, "Length2 element too small on primitive %d!\n", prim);
308 nb2 = 1;
309 } else if (nb2 < MinNbElementsOnLength) {
310 // Need to have at least MinNbElementsOnLength elements per wire primitive
311 nb2 = MinNbElementsOnLength;
312 double ellength = l2 / (double)nb2;
313 // which may not be possible if the length is small
314 if (ellength < MINDIST) {
315 // However, it is necessary to have at least one element!
316 nb2 = (int)(l2 / MINDIST);
317 if (nb2 < 1) {
318 fprintf(fMeshLog, "Length2 element very small on primitive %d!\n",
319 prim);
320 nb2 = 1;
321 }
322 }
323 } // else if nb2 < MinNbElementsOnLength
324
325 if (nb2 > MaxNbElementsOnLength) {
326 fprintf(fMeshLog, "Too many elements on Length2 of primitive %d!\n",
327 prim);
328 fprintf(fMeshLog, "Number of elements reduced to maximum allowed %d\n",
329 MaxNbElementsOnLength);
330 nb2 = MaxNbElementsOnLength;
331 }
332
333 *NbSegCoord1 = nb1;
334 *NbSegCoord2 = nb2;
335
336 fprintf(fMeshLog,
337 "Number of elements on surface primitive %d is %d X %d.\n\n", prim,
338 *NbSegCoord1, *NbSegCoord2);
339 } else { // number of discretization specified by the user
340 // Triangle primitives have their right angle on the vertex 1
341 double l1 = (XVertex[prim][0] - XVertex[prim][1]) *
342 (XVertex[prim][0] - XVertex[prim][1]) +
343 (YVertex[prim][0] - YVertex[prim][1]) *
344 (YVertex[prim][0] - YVertex[prim][1]) +
345 (ZVertex[prim][0] - ZVertex[prim][1]) *
346 (ZVertex[prim][0] - ZVertex[prim][1]);
347 l1 = sqrt(l1);
348 double l2 = (XVertex[prim][2] - XVertex[prim][1]) *
349 (XVertex[prim][2] - XVertex[prim][1]) +
350 (YVertex[prim][2] - YVertex[prim][1]) *
351 (YVertex[prim][2] - YVertex[prim][1]) +
352 (ZVertex[prim][2] - ZVertex[prim][1]) *
353 (ZVertex[prim][2] - ZVertex[prim][1]);
354 l2 = sqrt(l2);
355
356 if (l1 > l2) {
357 if (nb2 > nb1) {
358 // swap numbers to allow larger number to larger side
359 int tmpnb = nb1;
360 nb1 = nb2;
361 nb2 = tmpnb;
362 }
363 } // if l1 > l2
364
365 double ellength1 = l1 / (double)nb1;
366 double ellength2 = l2 / (double)nb2;
367
368 if (l1 < MINDIST) {
369 nb1 = 1;
370 fprintf(fMeshLog, "Fatal: Side length l1 too small! prim: %d\n", prim);
371 } else if (ellength1 < MINDIST) // element length more than twice MINDIST
372 {
373 nb1 = (int)(l1 / (2.0 * MINDIST));
374 if (nb1 < 1) {
375 nb1 = 1;
376 fprintf(fMeshLog, "Fatal: Side length l1 too small on primitive %d!\n",
377 prim);
378 }
379 } // if ellength1 < MINDIST
380
381 if (l2 < MINDIST) {
382 nb2 = 1;
383 fprintf(fMeshLog, "Fatal: Side length l2 too small! prim: %d\n", prim);
384 } else if (ellength2 < MINDIST) // element length more than twice MINDIST
385 {
386 nb2 = (int)(l2 / (2.0 * MINDIST));
387 if (nb2 < 1) {
388 nb2 = 1;
389 fprintf(fMeshLog, "Fatal: Side length l2 too small on primitive %d!\n",
390 prim);
391 }
392 } // if ellength2 < MINDIST
393
394 *NbSegCoord1 = nb1;
395 *NbSegCoord2 = nb2;
396
397 fprintf(fMeshLog,
398 "Number of elements on surface primitive %d is %d X %d.\n\n", prim,
399 *NbSegCoord1, *NbSegCoord2);
400 }
401
402 if ((nb1 > 0) && (nb2 > 0))
403 return 0;
404 else
405 return -1;
406} // AnalyzeSurface ends
MINDIST
#define MINDIST
Definition: Isles.h:17
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314
ZVertex
neBEMGLOBAL double ** ZVertex
Definition: neBEM.h:71
ElementLengthRqstd
neBEMGLOBAL double ElementLengthRqstd
Definition: neBEM.h:107
YVertex
neBEMGLOBAL double ** YVertex
Definition: neBEM.h:71
XVertex
neBEMGLOBAL double ** XVertex
Definition: neBEM.h:71
fMeshLog
neBEMGLOBAL FILE * fMeshLog
Definition: neBEM.h:114
MaxNbElementsOnLength
neBEMGLOBAL int MaxNbElementsOnLength
Definition: neBEM.h:105
MinNbElementsOnLength
neBEMGLOBAL int MinNbElementsOnLength
Definition: neBEM.h:103

Referenced by AnalyzePrimitive().

◆ AnalyzeWire()

int AnalyzeWire ( int  prim,
int *  NbSeg 
)

Definition at line 158 of file ReTriM.c.

158 {
159 int nb = *NbSeg;
160
161 if (nb < 1) // absurd! use the trio: target, min, max
162 {
163 double lWire = (XVertex[prim][1] - XVertex[prim][0]) *
164 (XVertex[prim][1] - XVertex[prim][0]) +
165 (YVertex[prim][1] - YVertex[prim][0]) *
166 (YVertex[prim][1] - YVertex[prim][0]) +
167 (ZVertex[prim][1] - ZVertex[prim][0]) *
168 (ZVertex[prim][1] - ZVertex[prim][0]);
169 lWire = sqrt(lWire);
170
171 nb = (int)(lWire / ElementLengthRqstd);
172
173 if ((nb > MinNbElementsOnLength) &&
174 (nb < MaxNbElementsOnLength)) { // nothing to be done
175 }
176 // Check whether the length of the wire primitive is long enough
177 else if (lWire < MINDIST) {
178 nb = 1;
179 fprintf(fMeshLog, "Wire element too small on primitive %d!\n", prim);
180 } // if lWire < MINDIST
181 // Need to have at least MinNbElementsOnLength elements per wire primitive
182 else if (nb < MinNbElementsOnLength) {
183 nb = MinNbElementsOnLength;
184 double ellength = lWire / (double)nb;
185 if (ellength <
186 MINDIST) // which may not be possible if the length is small
187 {
188 nb = (int)(lWire / MINDIST);
189 if (nb < 1) // However, it is necessary to have at least one element!
190 {
191 nb = 1;
192 fprintf(fMeshLog, "Wire element very small on primitive %d!\n", prim);
193 }
194 }
195 } // if nb < MinNbElementsOnLength
196 else if (nb > MaxNbElementsOnLength) {
197 nb = MaxNbElementsOnLength;
198 fprintf(fMeshLog, "Too many elements on wire primitive %d!\n", prim);
199 fprintf(fMeshLog, "Number of elements reduced to maximum allowed %d\n",
200 MaxNbElementsOnLength);
201 } // if nb > MaxNbElementsOnLength
202
203 *NbSeg = nb;
204
205 fprintf(fMeshLog, "Number of elements on wire primitive %d is %d.\n\n",
206 prim, *NbSeg);
207 } else { // number of dicretization specified by user
208 double lWire = (XVertex[prim][1] - XVertex[prim][0]) *
209 (XVertex[prim][1] - XVertex[prim][0]) +
210 (YVertex[prim][1] - YVertex[prim][0]) *
211 (YVertex[prim][1] - YVertex[prim][0]) +
212 (ZVertex[prim][1] - ZVertex[prim][0]) *
213 (ZVertex[prim][1] - ZVertex[prim][0]);
214 lWire = sqrt(lWire);
215
216 double ellength = lWire / (double)nb;
217
218 if (lWire < MINDIST) // at least one element is a necessity
219 {
220 nb = 1;
221 fprintf(fMeshLog, "Fatal: Wire element too small on primitive %d!\n",
222 prim);
223 } // if lWire < MINDIST
224 else if (ellength < MINDIST) // element length more than twice MINDIST
225 {
226 nb = (int)(lWire / (2.0 * MINDIST));
227 if (nb < 1) {
228 nb = 1;
229 fprintf(fMeshLog, "Fatal: Wire element too small on primitive %d!\n",
230 prim);
231 }
232 } // if ellength < MINDIST
233
234 *NbSeg = nb;
235
236 fprintf(fMeshLog, "Number of elements on wire primitive %d is %d.\n\n",
237 prim, *NbSeg);
238 }
239
240 if (nb)
241 return 0;
242 else
243 return -1;
244} // AnalyzeWire ends

Referenced by AnalyzePrimitive().

◆ BoundaryConditions()

int BoundaryConditions ( void  )

Definition at line 2121 of file ReTriM.c.

2121 {
2122 // Loop over all the elements of all the primitives.
2123 // Some of the primitives do not have a boundary condition to be satisfied
2124 // We'll omit these primitives but before doing that we need to know to
2125 // which primitive a particular element belongs to.
2126 // The RHS will also need to be modified to accommodate the presence of
2127 // floating conductors and charged (known) substances.
2128 // Looping on primitives (rather than elements) can save precious time!
2129 for (int ele = 1; ele <= NbElements; ++ele) {
2130 int prim = (EleArr + ele - 1)->PrimitiveNb;
2131 // Note this condition is pure geometry, not electric!
2132 switch (PrimType[prim]) {
2133 case 2:
2134 case 3: // for floating conductor and for dielectric-dielectric intfc
2135 case 4: // the BC is zero (may be modified due to known charges, later)
2136 (EleArr + ele - 1)->BC.Value = ApplPot[prim];
2137 break;
2138
2139 default:
2140 printf("Primitive out of range in BoundaryConditions ... returning\n");
2141 return (-1);
2142 } // switch on PrimType over
2143 } // ele loop over
2144
2145 return (0);
2146} // end of BoundaryConditions
EleArr
neBEMGLOBAL Element * EleArr
Definition: neBEM.h:169
ApplPot
neBEMGLOBAL double * ApplPot
Definition: neBEM.h:75
NbElements
neBEMGLOBAL int NbElements
Definition: neBEM.h:113
PrimType
neBEMGLOBAL int * PrimType
Definition: neBEM.h:67

Referenced by neBEMBoundaryConditions().

◆ DiscretizePolygon()

int DiscretizePolygon ( 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 at line 2106 of file ReTriM.c.

2111 {
2112 // Check inputs
2113 if ((NbSegX <= 0) || (NbSegZ <= 0)) {
2114 printf("segmentation input wrong in DiscretizePolygon ...\n");
2115 return -1;
2116 }
2117
2118 return (0);
2119} // end of DiscretizePolygon

◆ DiscretizeRectangle()

int DiscretizeRectangle ( 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 at line 1590 of file ReTriM.c.

1594 {
1595 int SurfParentObj, SurfEType;
1596 double SurfX, SurfY, SurfZ, SurfLX, SurfLZ;
1597 double SurfElX, SurfElY, SurfElZ, SurfElLX, SurfElLZ;
1598 DirnCosn3D PrimDirnCosn; // direction cosine of the current primitive
1599 char primstr[10];
1600 char gpElem[256], gpMesh[256];
1601 FILE *fPrim, *fElem, *fgpPrim, *fgpElem, *fgpMesh;
1602
1603 // Check inputs
1604 if ((NbSegX <= 0) || (NbSegZ <= 0)) {
1605 printf("segmentation input wrong in DiscretizeRectangle ...\n");
1606 return -1;
1607 }
1608
1609 if (OptPrintVertexAndNormal) {
1610 printf("nvertex: %d\n", nvertex);
1611 for (int vert = 0; vert < nvertex; ++vert) {
1612 printf("vert: %d, x: %lg, y: %lg, z: %lg\n", vert, xvert[vert],
1613 yvert[vert], zvert[vert]);
1614 }
1615 printf("Normal: %lg, %lg, %lg\n", xnorm, ynorm, znorm);
1616 } // if OptPrintVertexAndNormal
1617
1618 // necessary for separating filenames
1619 sprintf(primstr, "%d", prim);
1620
1621 // in order to avoid warning messages
1622 fPrim = NULL;
1623 fElem = NULL;
1624 fgpMesh = NULL;
1625 fgpElem = NULL;
1626
1627 // Get volume information, to begin with
1628 // int shape, material, boundarytype;
1629 // double eps, potential, charge;
1630 // neBEMVolumeDescription(volref1, &shape, &material,
1631 // &eps, &potential, &charge, &boundarytype);
1632
1633 // compute all the properties of this surface
1634 SurfParentObj = 1;
1635 SurfEType = inttype;
1636 if (SurfEType == 0) {
1637 printf("Wrong SurfEType for prim %d\n", prim);
1638 exit(-1);
1639 }
1640 // centroid of the local coordinate system
1641 SurfX = 0.25 * (xvert[0] + xvert[1] + xvert[2] + xvert[3]);
1642 SurfY = 0.25 * (yvert[0] + yvert[1] + yvert[2] + yvert[3]);
1643 SurfZ = 0.25 * (zvert[0] + zvert[1] + zvert[2] + zvert[3]);
1644 // lengths of the sides
1645 SurfLX = sqrt((xvert[1] - xvert[0]) * (xvert[1] - xvert[0]) +
1646 (yvert[1] - yvert[0]) * (yvert[1] - yvert[0]) +
1647 (zvert[1] - zvert[0]) * (zvert[1] - zvert[0]));
1648 SurfLZ = sqrt((xvert[2] - xvert[1]) * (xvert[2] - xvert[1]) +
1649 (yvert[2] - yvert[1]) * (yvert[2] - yvert[1]) +
1650 (zvert[2] - zvert[1]) * (zvert[2] - zvert[1]));
1651 // Direction cosines
1652 PrimDirnCosn.XUnit.X = (xvert[1] - xvert[0]) / SurfLX;
1653 PrimDirnCosn.XUnit.Y = (yvert[1] - yvert[0]) / SurfLX;
1654 PrimDirnCosn.XUnit.Z = (zvert[1] - zvert[0]) / SurfLX;
1655 PrimDirnCosn.YUnit.X = xnorm;
1656 PrimDirnCosn.YUnit.Y = ynorm;
1657 PrimDirnCosn.YUnit.Z = znorm;
1658 PrimDirnCosn.ZUnit =
1659 Vector3DCrossProduct(&PrimDirnCosn.XUnit, &PrimDirnCosn.YUnit);
1660
1661 // primitive direction cosine assignments
1662 PrimDC[prim].XUnit.X = PrimDirnCosn.XUnit.X;
1663 PrimDC[prim].XUnit.Y = PrimDirnCosn.XUnit.Y;
1664 PrimDC[prim].XUnit.Z = PrimDirnCosn.XUnit.Z;
1665 PrimDC[prim].YUnit.X = PrimDirnCosn.YUnit.X;
1666 PrimDC[prim].YUnit.Y = PrimDirnCosn.YUnit.Y;
1667 PrimDC[prim].YUnit.Z = PrimDirnCosn.YUnit.Z;
1668 PrimDC[prim].ZUnit.X = PrimDirnCosn.ZUnit.X;
1669 PrimDC[prim].ZUnit.Y = PrimDirnCosn.ZUnit.Y;
1670 PrimDC[prim].ZUnit.Z = PrimDirnCosn.ZUnit.Z;
1671
1672 // primitive origin: also the barcenter for a rectangular element
1673 PrimOriginX[prim] = SurfX;
1674 PrimOriginY[prim] = SurfY;
1675 PrimOriginZ[prim] = SurfZ;
1676 PrimLX[prim] = SurfLX;
1677 PrimLZ[prim] = SurfLZ;
1678
1679 double SurfLambda = lambda;
1680 double SurfV = potential;
1681
1682 // file output for a primitive
1683 if (OptPrimitiveFiles) {
1684 char OutPrim[256];
1685 strcpy(OutPrim, ModelOutDir);
1686 strcat(OutPrim, "/Primitives/Primitive");
1687 strcat(OutPrim, primstr);
1688 strcat(OutPrim, ".out");
1689 fPrim = fopen(OutPrim, "w");
1690 if (fPrim == NULL) {
1691 neBEMMessage("DiscretizeRectangle - OutPrim");
1692 return -1;
1693 }
1694 fprintf(fPrim, "#prim: %d, nvertex: %d\n", prim, nvertex);
1695 fprintf(fPrim, "Node1: %lg\t%lg\t%lg\n", xvert[0], yvert[0], zvert[0]);
1696 fprintf(fPrim, "Node2: %lg\t%lg\t%lg\n", xvert[1], yvert[1], zvert[1]);
1697 fprintf(fPrim, "Node3: %lg\t%lg\t%lg\n", xvert[2], yvert[2], zvert[2]);
1698 fprintf(fPrim, "Node4: %lg\t%lg\t%lg\n", xvert[3], yvert[3], zvert[3]);
1699 fprintf(fPrim, "PrimOrigin: %lg\t%lg\t%lg\n", PrimOriginX[prim],
1700 PrimOriginY[prim], PrimOriginZ[prim]);
1701 fprintf(fPrim, "Primitive lengths: %lg\t%lg\n", PrimLX[prim], PrimLZ[prim]);
1702 fprintf(fPrim, "Norm: %lg\t%lg\t%lg\n", xnorm, ynorm, znorm);
1703 fprintf(fPrim, "#volref1: %d, volref2: %d\n", volref1, volref2);
1704 fprintf(fPrim, "#NbSegX: %d, NbSegZ: %d\n", NbSegX, NbSegZ);
1705 fprintf(fPrim, "#ParentObj: %d\tEType: %d\n", SurfParentObj, SurfEType);
1706 fprintf(fPrim, "#SurfX\tSurfY\tSurfZ\tSurfLZ\tSurfLZ\n");
1707 fprintf(fPrim, "%lg\t%lg\t%lg\t%lg\t%lg\n", SurfX, SurfY, SurfZ, SurfLX,
1708 SurfLZ);
1709 // fprintf(fPrim, "#SurfRX: %lg\tSurfRY: %lg\tSurfRZ: %lg\n",
1710 // SurfRX, SurfRY, SurfRZ);
1711 fprintf(fPrim, "#DirnCosn: \n");
1712 fprintf(fPrim, "%lg, %lg, %lg\n", PrimDirnCosn.XUnit.X,
1713 PrimDirnCosn.XUnit.Y, PrimDirnCosn.XUnit.Z);
1714 fprintf(fPrim, "%lg, %lg, %lg\n", PrimDirnCosn.YUnit.X,
1715 PrimDirnCosn.YUnit.Y, PrimDirnCosn.YUnit.Z);
1716 fprintf(fPrim, "%lg, %lg, %lg\n", PrimDirnCosn.ZUnit.X,
1717 PrimDirnCosn.ZUnit.Y, PrimDirnCosn.ZUnit.Z);
1718 fprintf(fPrim, "#SurfLambda: %lg\tSurfV: %lg\n", SurfLambda, SurfV);
1719 } // if OptPrimitiveFiles
1720
1721 // necessary for gnuplot
1722 if (OptGnuplot && OptGnuplotPrimitives) {
1723 char gpPrim[256];
1724 strcpy(gpPrim, MeshOutDir);
1725 strcat(gpPrim, "/GViewDir/gpPrim");
1726 strcat(gpPrim, primstr);
1727 strcat(gpPrim, ".out");
1728 fgpPrim = fopen(gpPrim, "w");
1729 if (fgpPrim == NULL) {
1730 neBEMMessage("DiscretizeRectangle - OutgpPrim");
1731 return -1;
1732 }
1733 fprintf(fgpPrim, "%g\t%g\t%g\n\n", xvert[0], yvert[0], zvert[0]);
1734 fprintf(fgpPrim, "%g\t%g\t%g\n\n", xvert[1], yvert[1], zvert[1]);
1735 fprintf(fgpPrim, "%g\t%g\t%g\n\n", xvert[2], yvert[2], zvert[2]);
1736 fprintf(fgpPrim, "%g\t%g\t%g\n\n", xvert[3], yvert[3], zvert[3]);
1737 fprintf(fgpPrim, "%g\t%g\t%g\n", xvert[0], yvert[0], zvert[0]);
1738 fclose(fgpPrim);
1739
1740 if (prim == 1)
1741 fprintf(fgnuPrim, " '%s\' w l", gpPrim);
1742 else
1743 fprintf(fgnuPrim, ", \\\n \'%s\' w l", gpPrim);
1744 } // if OptGnuplot && OptGnuplotPrimitives
1745
1746 // file outputs for elements on primitive
1747 if (OptElementFiles) {
1748 char OutElem[256];
1749 strcpy(OutElem, MeshOutDir);
1750 strcat(OutElem, "/Elements/ElemOnPrim");
1751 strcat(OutElem, primstr);
1752 strcat(OutElem, ".out");
1753 fElem = fopen(OutElem, "w");
1754 // assert(fElem != NULL);
1755 if (fElem == NULL) {
1756 neBEMMessage("DiscretizeRectangle - OutElem");
1757 return -1;
1758 }
1759 } // if OptElementFiles
1760
1761 // gnuplot friendly file outputs for elements on primitive
1762 if (OptGnuplot && OptGnuplotElements) {
1763 strcpy(gpElem, MeshOutDir);
1764 strcat(gpElem, "/GViewDir/gpElemOnPrim");
1765 strcat(gpElem, primstr);
1766 strcat(gpElem, ".out");
1767 fgpElem = fopen(gpElem, "w");
1768 // assert(fgpElem != NULL);
1769 if (fgpElem == NULL) {
1770 neBEMMessage("DiscretizeRectangle - OutgpElem");
1771 if (fElem) fclose(fElem);
1772 return -1;
1773 }
1774 // gnuplot friendly file outputs for elements on primitive
1775 strcpy(gpMesh, MeshOutDir);
1776 strcat(gpMesh, "/GViewDir/gpMeshOnPrim");
1777 strcat(gpMesh, primstr);
1778 strcat(gpMesh, ".out");
1779 fgpMesh = fopen(gpMesh, "w");
1780 if (fgpMesh == NULL) {
1781 neBEMMessage("DiscretizeRectangle - OutgpMesh");
1782 fclose(fgpElem);
1783 return -1;
1784 }
1785 } // if OptGnuplot && OptElements
1786
1787 // Compute element positions (CGs) in the primitive local coordinate system.
1788 // Then map these CGs to the global coordinate system.
1789 // (xav, 0, zav) is the CG of an element wrt the primitive coordinate system
1790 // From this, we find the offset of the element from the primitive CG in the
1791 // global coordinate system by just rotating the displacement vector
1792 // (0,0,0) -> (xav,0,zav) using the known DCM of the surface.
1793 // This rotated vector (now in the global system) is then added to the
1794 // position vector of the primitive CG (always in the global system) to get
1795 // the CG of the element in the global system. Make sure that the larger
1796 // number is being used to discretize the longer side - can be OverSmart in
1797 // certain cases - there may be cases where smaller number of elements suffice
1798 // for a longer side
1799 if (NbSegX == NbSegZ) {
1800 SurfElLX = SurfLX / NbSegX; // element sizes
1801 SurfElLZ = SurfLZ / NbSegZ;
1802 } else if (NbSegX > NbSegZ) {
1803 if (SurfLX > SurfLZ) {
1804 SurfElLX = SurfLX / NbSegX; // element sizes
1805 SurfElLZ = SurfLZ / NbSegZ;
1806 } else // interchange NbSegX and NbSegZ
1807 {
1808 int tmp = NbSegZ;
1809 NbSegZ = NbSegX;
1810 NbSegX = tmp;
1811 SurfElLX = SurfLX / NbSegX; // element sizes
1812 SurfElLZ = SurfLZ / NbSegZ;
1813 }
1814 } // NbSegX > NbSegZ
1815 else // NbSegX < NbSegZ
1816 {
1817 if (SurfLX < SurfLZ) {
1818 SurfElLX = SurfLX / NbSegX; // element sizes
1819 SurfElLZ = SurfLZ / NbSegZ;
1820 } else // interchange NbSegX and NbSegZ
1821 {
1822 int tmp = NbSegZ;
1823 NbSegZ = NbSegX;
1824 NbSegX = tmp;
1825 SurfElLX = SurfLX / NbSegX; // element sizes
1826 SurfElLZ = SurfLZ / NbSegZ;
1827 }
1828 }
1829
1830 // Analysis of element aspect ratio.
1831 // Note that we can afford only to reduce the total number of elements.
1832 // Otherwise, we'll have to realloc `EleArr' array.
1833 // Later, we'll make the corrections such that the total number of elements
1834 // remain close to the originally intended value.
1835 double AR = SurfElLX / SurfElLZ; // indicative element aspect ratio
1836 if (OptPrimitiveFiles) {
1837 fprintf(fPrim,
1838 "Using the input, the aspect ratio of the elements on prim: %d\n",
1839 prim);
1840 fprintf(fPrim,
1841 "NbSegX: %d, SurfElLX: %lg, NbSegZ: %d, SurfElLZ: %lg, AR: %lg\n",
1842 NbSegX, SurfElLX, NbSegZ, SurfElLZ, AR);
1843 }
1844 if (AR > 10.0) // NbSegZ needs to be reduced
1845 {
1846 double tmpElLZ = SurfElLX / 10.0;
1847 NbSegZ = (int)(SurfLZ / tmpElLZ);
1848 if (NbSegZ <= 0) NbSegZ = 1;
1849 SurfElLZ = SurfLZ / NbSegZ;
1850 }
1851 if (AR < 0.1) // NbSegX need to be reduced
1852 {
1853 double tmpElLX = SurfElLZ * 0.1;
1854 NbSegX = (int)(SurfLX / tmpElLX);
1855 if (NbSegX <= 0) NbSegX = 1;
1856 SurfElLX = SurfLX / NbSegX;
1857 }
1858 AR = SurfElLX / SurfElLZ; // indicative element aspect ratio
1859 if (OptPrimitiveFiles) {
1860 fprintf(fPrim, "After analyzing the aspect ratio of the elements\n");
1861 fprintf(fPrim,
1862 "NbSegX: %d, SurfElLX: %lg, NbSegZ: %d, SurfElLZ: %lg, AR: %lg\n",
1863 NbSegX, SurfElLX, NbSegZ, SurfElLZ, AR);
1864 }
1865
1866 double x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3, xav, zav;
1867 ElementBgn[prim] = EleCntr + 1;
1868 for (int i = 1; i <= NbSegX; ++i) {
1869 x1 = -SurfLX / 2.0 +
1870 (double)(i - 1) * SurfElLX; // assuming centroid at 0,0,0
1871 x2 = -SurfLX / 2.0 + (double)(i)*SurfElLX;
1872 xav = 0.5 * (x1 + x2);
1873
1874 for (int k = 1; k <= NbSegZ; ++k) {
1875 z1 = -SurfLZ / 2.0 + (double)(k - 1) * SurfElLZ;
1876 z2 = -SurfLZ / 2.0 + (double)(k)*SurfElLZ;
1877 zav = 0.5 * (z1 + z2);
1878
1879 { // Separate block for position rotation - local2global
1880 Point3D localDisp, globalDisp;
1881
1882 localDisp.X = xav;
1883 localDisp.Y = 0.0;
1884 localDisp.Z = zav;
1885 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1886 SurfElX = SurfX + globalDisp.X;
1887 SurfElY = SurfY + globalDisp.Y;
1888 SurfElZ = SurfZ + globalDisp.Z;
1889 } // vector rotation over
1890
1891 // Assign element values and write in the file
1892 ++EleCntr;
1893 if (EleCntr > NbElements) {
1894 neBEMMessage("DiscretizeRectangle - EleCntr more than NbElements!");
1895 if (fgpMesh) fclose(fgpMesh);
1896 return -1;
1897 }
1898
1899 (EleArr + EleCntr - 1)->DeviceNb =
1900 1; // At present, there is only one device
1901 (EleArr + EleCntr - 1)->ComponentNb = SurfParentObj;
1902 (EleArr + EleCntr - 1)->PrimitiveNb = prim;
1903 (EleArr + EleCntr - 1)->Id = EleCntr;
1904 (EleArr + EleCntr - 1)->G.Type = 4; // rectangular here
1905 (EleArr + EleCntr - 1)->G.Origin.X = SurfElX;
1906 (EleArr + EleCntr - 1)->G.Origin.Y = SurfElY;
1907 (EleArr + EleCntr - 1)->G.Origin.Z = SurfElZ;
1908 (EleArr + EleCntr - 1)->G.LX = SurfElLX;
1909 (EleArr + EleCntr - 1)->G.LZ = SurfElLZ;
1910 (EleArr + EleCntr - 1)->G.dA =
1911 (EleArr + EleCntr - 1)->G.LX * (EleArr + EleCntr - 1)->G.LZ;
1912 (EleArr + EleCntr - 1)->G.DC.XUnit.X = PrimDirnCosn.XUnit.X;
1913 (EleArr + EleCntr - 1)->G.DC.XUnit.Y = PrimDirnCosn.XUnit.Y;
1914 (EleArr + EleCntr - 1)->G.DC.XUnit.Z = PrimDirnCosn.XUnit.Z;
1915 (EleArr + EleCntr - 1)->G.DC.YUnit.X = PrimDirnCosn.YUnit.X;
1916 (EleArr + EleCntr - 1)->G.DC.YUnit.Y = PrimDirnCosn.YUnit.Y;
1917 (EleArr + EleCntr - 1)->G.DC.YUnit.Z = PrimDirnCosn.YUnit.Z;
1918 (EleArr + EleCntr - 1)->G.DC.ZUnit.X = PrimDirnCosn.ZUnit.X;
1919 (EleArr + EleCntr - 1)->G.DC.ZUnit.Y = PrimDirnCosn.ZUnit.Y;
1920 (EleArr + EleCntr - 1)->G.DC.ZUnit.Z = PrimDirnCosn.ZUnit.Z;
1921 (EleArr + EleCntr - 1)->E.Type = SurfEType;
1922 (EleArr + EleCntr - 1)->E.Lambda = SurfLambda;
1923 (EleArr + EleCntr - 1)->Solution = 0.0;
1924 (EleArr + EleCntr - 1)->Assigned = charge;
1925 (EleArr + EleCntr - 1)->BC.NbOfBCs = 1; // assume one BC per element
1926 (EleArr + EleCntr - 1)->BC.CollPt.X = (EleArr + EleCntr - 1)->G.Origin.X;
1927 (EleArr + EleCntr - 1)->BC.CollPt.Y = (EleArr + EleCntr - 1)->G.Origin.Y;
1928 (EleArr + EleCntr - 1)->BC.CollPt.Z = (EleArr + EleCntr - 1)->G.Origin.Z;
1929
1930 // find element vertices
1931 // 1) displacement vector in the ECS is first identified
1932 // 2) this vector is transformed to the GCS
1933 // 3) the global displacement vector, when added to the centroid in GCS,
1934 // gives the node positions in GCS.
1935 x0 = -0.5 *
1936 (EleArr + EleCntr - 1)->G.LX; // xyz displacement of node wrt ECS
1937 y0 = 0.0;
1938 z0 = -0.5 * (EleArr + EleCntr - 1)->G.LZ;
1939 { // Separate block for position rotation - local2global
1940 Point3D localDisp, globalDisp;
1941
1942 localDisp.X = x0;
1943 localDisp.Y = y0;
1944 localDisp.Z = z0; // displacement in GCS
1945 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1946 x0 = (EleArr + EleCntr - 1)->G.Origin.X +
1947 globalDisp.X; // xyz position in GCS
1948 y0 = (EleArr + EleCntr - 1)->G.Origin.Y + globalDisp.Y;
1949 z0 = (EleArr + EleCntr - 1)->G.Origin.Z + globalDisp.Z;
1950 } // position rotation over
1951
1952 x1 = 0.5 * (EleArr + EleCntr - 1)->G.LX;
1953 y1 = 0.0;
1954 z1 = -0.5 * (EleArr + EleCntr - 1)->G.LZ;
1955 { // Separate block for position rotation - local2global
1956 Point3D localDisp, globalDisp;
1957
1958 localDisp.X = x1;
1959 localDisp.Y = y1;
1960 localDisp.Z = z1;
1961 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1962 x1 = (EleArr + EleCntr - 1)->G.Origin.X + globalDisp.X;
1963 y1 = (EleArr + EleCntr - 1)->G.Origin.Y + globalDisp.Y;
1964 z1 = (EleArr + EleCntr - 1)->G.Origin.Z + globalDisp.Z;
1965 } // position rotation over
1966
1967 x2 = 0.5 * (EleArr + EleCntr - 1)->G.LX;
1968 y2 = 0.0;
1969 z2 = 0.5 * (EleArr + EleCntr - 1)->G.LZ;
1970 { // Separate block for position rotation - local2global
1971 Point3D localDisp, globalDisp;
1972
1973 localDisp.X = x2;
1974 localDisp.Y = y2;
1975 localDisp.Z = z2;
1976 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1977 x2 = (EleArr + EleCntr - 1)->G.Origin.X + globalDisp.X;
1978 y2 = (EleArr + EleCntr - 1)->G.Origin.Y + globalDisp.Y;
1979 z2 = (EleArr + EleCntr - 1)->G.Origin.Z + globalDisp.Z;
1980 } // position rotation over
1981
1982 x3 = -0.5 * (EleArr + EleCntr - 1)->G.LX;
1983 y3 = 0.0;
1984 z3 = 0.5 * (EleArr + EleCntr - 1)->G.LZ;
1985 { // Separate block for position rotation - local2global
1986 Point3D localDisp, globalDisp;
1987
1988 localDisp.X = x3;
1989 localDisp.Y = y3;
1990 localDisp.Z = z3;
1991 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1992 x3 = (EleArr + EleCntr - 1)->G.Origin.X + globalDisp.X;
1993 y3 = (EleArr + EleCntr - 1)->G.Origin.Y + globalDisp.Y;
1994 z3 = (EleArr + EleCntr - 1)->G.Origin.Z + globalDisp.Z;
1995 } // position rotation over
1996
1997 // assign vertices of the element
1998 (EleArr + EleCntr - 1)->G.Vertex[0].X = x0;
1999 (EleArr + EleCntr - 1)->G.Vertex[0].Y = y0;
2000 (EleArr + EleCntr - 1)->G.Vertex[0].Z = z0;
2001 (EleArr + EleCntr - 1)->G.Vertex[1].X = x1;
2002 (EleArr + EleCntr - 1)->G.Vertex[1].Y = y1;
2003 (EleArr + EleCntr - 1)->G.Vertex[1].Z = z1;
2004 (EleArr + EleCntr - 1)->G.Vertex[2].X = x2;
2005 (EleArr + EleCntr - 1)->G.Vertex[2].Y = y2;
2006 (EleArr + EleCntr - 1)->G.Vertex[2].Z = z2;
2007 (EleArr + EleCntr - 1)->G.Vertex[3].X = x3;
2008 (EleArr + EleCntr - 1)->G.Vertex[3].Y = y3;
2009 (EleArr + EleCntr - 1)->G.Vertex[3].Z = z3;
2010
2011 if (OptElementFiles) {
2012 fprintf(fElem, "##Element Counter: %d\n", EleCntr);
2013 fprintf(fElem, "#DevNb\tCompNb\tPrimNb\tId\n");
2014 fprintf(fElem, "%d\t%d\t%d\t%d\n", (EleArr + EleCntr - 1)->DeviceNb,
2015 (EleArr + EleCntr - 1)->ComponentNb,
2016 (EleArr + EleCntr - 1)->PrimitiveNb,
2017 (EleArr + EleCntr - 1)->Id);
2018 fprintf(fElem, "#GType\tX\tY\tZ\tLX\tLZ\tdA\n");
2019 fprintf(
2020 fElem, "%d\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\n",
2021 (EleArr + EleCntr - 1)->G.Type, (EleArr + EleCntr - 1)->G.Origin.X,
2022 (EleArr + EleCntr - 1)->G.Origin.Y,
2023 (EleArr + EleCntr - 1)->G.Origin.Z, (EleArr + EleCntr - 1)->G.LX,
2024 (EleArr + EleCntr - 1)->G.LZ, (EleArr + EleCntr - 1)->G.dA);
2025 fprintf(fElem, "#DirnCosn: \n");
2026 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.XUnit.X,
2027 (EleArr + EleCntr - 1)->G.DC.XUnit.Y,
2028 (EleArr + EleCntr - 1)->G.DC.XUnit.Z);
2029 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.YUnit.X,
2030 (EleArr + EleCntr - 1)->G.DC.YUnit.Y,
2031 (EleArr + EleCntr - 1)->G.DC.YUnit.Z);
2032 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.ZUnit.X,
2033 (EleArr + EleCntr - 1)->G.DC.ZUnit.Y,
2034 (EleArr + EleCntr - 1)->G.DC.ZUnit.Z);
2035 fprintf(fElem, "#EType\tLambda\n");
2036 fprintf(fElem, "%d\t%lg\n", (EleArr + EleCntr - 1)->E.Type,
2037 (EleArr + EleCntr - 1)->E.Lambda);
2038 fprintf(fElem, "#NbBCs\tCPX\tCPY\tCPZ\tValue\n");
2039 fprintf(fElem, "%d\t%lg\t%lg\t%lg\t%lg\n",
2040 (EleArr + EleCntr - 1)->BC.NbOfBCs,
2041 (EleArr + EleCntr - 1)->BC.CollPt.X,
2042 (EleArr + EleCntr - 1)->BC.CollPt.Y,
2043 (EleArr + EleCntr - 1)->BC.CollPt.Z,
2044 (EleArr + EleCntr - 1)->BC.Value);
2045 } // if OptElementFiles
2046
2047 // mark centroid
2048 if (OptGnuplot && OptGnuplotElements) {
2049 fprintf(fgpElem, "%g\t%g\t%g\n", (EleArr + EleCntr - 1)->BC.CollPt.X,
2050 (EleArr + EleCntr - 1)->BC.CollPt.Y,
2051 (EleArr + EleCntr - 1)->BC.CollPt.Z);
2052 } // if OptGnuplot && OptGnuplotElements
2053
2054 if (OptGnuplot && OptGnuplotElements) {
2055 fprintf(fgpMesh, "%g\t%g\t%g\n", x0, y0, z0);
2056 fprintf(fgpMesh, "%g\t%g\t%g\n", x1, y1, z1);
2057 fprintf(fgpMesh, "%g\t%g\t%g\n", x2, y2, z2);
2058 fprintf(fgpMesh, "%g\t%g\t%g\n", x3, y3, z3);
2059 fprintf(fgpMesh, "%g\t%g\t%g\n\n", x0, y0, z0);
2060 /*
2061 fprintf(fgpMesh, "%g\t%g\t%g\t%d\n", x0, y0, z0, 1);
2062 fprintf(fgpMesh, "%g\t%g\t%g\t%d\n\n", x1, y1, z1, 2);
2063 fprintf(fgpMesh, "%g\t%g\t%g\t%d\n", x2, y2, z2, 1);
2064 fprintf(fgpMesh, "%g\t%g\t%g\t%d\n\n\n", x2, y2, z2, 2);
2065 fprintf(fgpMesh, "%g\t%g\t%g\t%d\n", x0, y0, z0, 1);
2066 fprintf(fgpMesh, "%g\t%g\t%g\t%d\n\n", x3, y3, z3, 2);
2067 fprintf(fgpMesh, "%g\t%g\t%g\t%d\n", x2, y2, z2, 1);
2068 fprintf(fgpMesh, "%g\t%g\t%g\t%d\n\n\n", x2, y2, z2, 2);
2069 */
2070 } // if(OptGnuplot && OptGnuplotElements)
2071 } // for k
2072 } // for i
2073 ElementEnd[prim] = EleCntr;
2074 NbElmntsOnPrim[prim] = ElementEnd[prim] - ElementBgn[prim] + 1;
2075
2076 if (OptPrimitiveFiles) {
2077 fprintf(fPrim, "Element begin: %d, Element end: %d\n", ElementBgn[prim],
2078 ElementEnd[prim]);
2079 fprintf(fPrim, "Number of elements on primitive: %d\n",
2080 NbElmntsOnPrim[prim]);
2081 fclose(fPrim);
2082 }
2083
2084 if (OptElementFiles) fclose(fElem);
2085
2086 if (OptGnuplot && OptGnuplotElements) {
2087 if (prim == 1)
2088 fprintf(fgnuElem, " '%s\' w p", gpElem);
2089 else
2090 fprintf(fgnuElem, ", \\\n \'%s\' w p", gpElem);
2091 if (prim == 1) {
2092 fprintf(fgnuMesh, " '%s\' w l", gpMesh);
2093 fprintf(fgnuMesh, ", \\\n \'%s\' w p ps 1", gpElem);
2094 } else {
2095 fprintf(fgnuMesh, ", \\\n \'%s\' w l", gpMesh);
2096 fprintf(fgnuMesh, ", \\\n \'%s\' w p ps 1", gpElem);
2097 }
2098
2099 fclose(fgpElem);
2100 fclose(fgpMesh);
2101 }
2102
2103 return (0);
2104} // end of DiscretizeRectangle
RotatePoint3D
Point3D RotatePoint3D(Point3D *A, DirnCosn3D *DC, int Sense)
Definition: Vector.c:339
Vector3DCrossProduct
Vector3D Vector3DCrossProduct(Vector3D *A, Vector3D *B)
Definition: Vector.c:246
local2global
#define local2global
Definition: Vector.h:14
OptPrimitiveFiles
INTFACEGLOBAL int OptPrimitiveFiles
Definition: neBEMInterface.h:62
fgnuPrim
INTFACEGLOBAL FILE * fgnuPrim
Definition: neBEMInterface.h:55
fgnuElem
INTFACEGLOBAL FILE * fgnuElem
Definition: neBEMInterface.h:55
fgnuMesh
INTFACEGLOBAL FILE * fgnuMesh
Definition: neBEMInterface.h:55
OptGnuplotPrimitives
INTFACEGLOBAL int OptGnuplotPrimitives
Definition: neBEMInterface.h:54
OptPrintVertexAndNormal
INTFACEGLOBAL int OptPrintVertexAndNormal
Definition: neBEMInterface.h:62
OptElementFiles
INTFACEGLOBAL int OptElementFiles
Definition: neBEMInterface.h:98
OptGnuplot
INTFACEGLOBAL int OptGnuplot
Definition: neBEMInterface.h:54
OptGnuplotElements
INTFACEGLOBAL int OptGnuplotElements
Definition: neBEMInterface.h:54
NbElmntsOnPrim
neBEMGLOBAL int * NbElmntsOnPrim
Definition: neBEM.h:94
EleCntr
neBEMGLOBAL int EleCntr
Definition: neBEM.h:112
MeshOutDir
neBEMGLOBAL char MeshOutDir[256]
Definition: neBEM.h:263
PrimOriginZ
neBEMGLOBAL double * PrimOriginZ
Definition: neBEM.h:73
PrimOriginY
neBEMGLOBAL double * PrimOriginY
Definition: neBEM.h:73
PrimDC
neBEMGLOBAL DirnCosn3D * PrimDC
Definition: neBEM.h:74
ModelOutDir
neBEMGLOBAL char ModelOutDir[256]
Definition: neBEM.h:262
ElementEnd
neBEMGLOBAL int * ElementEnd
Definition: neBEM.h:94
PrimLX
neBEMGLOBAL double * PrimLX
Definition: neBEM.h:72
Solution
neBEMGLOBAL double * Solution
Definition: neBEM.h:237
ElementBgn
neBEMGLOBAL int * ElementBgn
Definition: neBEM.h:94
PrimOriginX
neBEMGLOBAL double * PrimOriginX
Definition: neBEM.h:73
PrimLZ
neBEMGLOBAL double * PrimLZ
Definition: neBEM.h:72
DirnCosn3D::ZUnit
Vector3D ZUnit
Definition: Vector.h:38
DirnCosn3D::YUnit
Vector3D YUnit
Definition: Vector.h:37
DirnCosn3D::XUnit
Vector3D XUnit
Definition: Vector.h:36
Point3D
Definition: Vector.h:21
Point3D::X
double X
Definition: Vector.h:22
Point3D::Z
double Z
Definition: Vector.h:24
Point3D::Y
double Y
Definition: Vector.h:23
Vector3D::Z
double Z
Definition: Vector.h:31
Vector3D::Y
double Y
Definition: Vector.h:30
Vector3D::X
double X
Definition: Vector.h:29

Referenced by SurfaceElements().

◆ DiscretizeTriangle()

int DiscretizeTriangle ( 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 at line 788 of file ReTriM.c.

791 {
792 int SurfParentObj, SurfEType;
793 double SurfX, SurfY, SurfZ, SurfLX, SurfLZ;
794 double SurfElX, SurfElY, SurfElZ, SurfElLX, SurfElLZ;
795 DirnCosn3D PrimDirnCosn; // direction cosine of the current primitive
796 char primstr[10];
797 char gpElem[256], gpMesh[256];
798 FILE *fPrim, *fElem, *fgpPrim, *fgpElem, *fgpMesh;
799
800 // Check inputs
801 if ((NbSegX <= 0) || (NbSegZ <= 0)) {
802 printf("segmentation input wrong in DiscretizeTriangle ...\n");
803 return -1;
804 }
805
806 if (OptPrintVertexAndNormal) {
807 printf("nvertex: %d\n", nvertex);
808 for (int vert = 0; vert < nvertex; ++vert) {
809 printf("vert: %d, x: %lg, y: %lg, z: %lg\n", vert, xvert[vert],
810 yvert[vert], zvert[vert]);
811 }
812 printf("Normal: %lg, %lg, %lg\n", xnorm, ynorm, znorm);
813 } // if OptPrintVertexAndNormal
814
815 // necessary for separating filenames
816 sprintf(primstr, "%d", prim);
817
818 // in order to avoid warning messages
819 fPrim = NULL;
820 fElem = NULL;
821 fgpMesh = NULL;
822 fgpElem = NULL;
823
824 // Compute all the properties of this surface
825 // Boundary types from 1 to 7 have been defined
826 SurfParentObj = 1;
827 SurfEType = inttype;
828 if ((SurfEType <= 0) || (SurfEType >= 8)) {
829 printf("Wrong SurfEType for prim %d\n", prim);
830 exit(-1);
831 }
832 // Origin of the local coordinate center is at the right angle corner
833 SurfX = xvert[1];
834 SurfY = yvert[1];
835 SurfZ = zvert[1];
836
837 // Find the proper direction cosines first - little more tricky that in the
838 // rectangular case
839 int flagDC = 0; // DC has not been ascertained as yet
840 // We begin with trial 1: one of the possible orientations
841 // Intially, the lengths are necessary
842 // lengths of the sides - note that right angle corner is [1]-th element
843 SurfLX = sqrt((xvert[0] - xvert[1]) * (xvert[0] - xvert[1]) +
844 (yvert[0] - yvert[1]) * (yvert[0] - yvert[1]) +
845 (zvert[0] - zvert[1]) * (zvert[0] - zvert[1]));
846 SurfLZ = sqrt((xvert[2] - xvert[1]) * (xvert[2] - xvert[1]) +
847 (yvert[2] - yvert[1]) * (yvert[2] - yvert[1]) +
848 (zvert[2] - zvert[1]) * (zvert[2] - zvert[1]));
849 // Direction cosines - note that right angle corner is [1]-th element
850 PrimDirnCosn.XUnit.X = (xvert[0] - xvert[1]) / SurfLX;
851 PrimDirnCosn.XUnit.Y = (yvert[0] - yvert[1]) / SurfLX;
852 PrimDirnCosn.XUnit.Z = (zvert[0] - zvert[1]) / SurfLX;
853 PrimDirnCosn.ZUnit.X = (xvert[2] - xvert[1]) / SurfLZ;
854 PrimDirnCosn.ZUnit.Y = (yvert[2] - yvert[1]) / SurfLZ;
855 PrimDirnCosn.ZUnit.Z = (zvert[2] - zvert[1]) / SurfLZ;
856 PrimDirnCosn.YUnit =
857 Vector3DCrossProduct(&PrimDirnCosn.ZUnit, &PrimDirnCosn.XUnit);
858 if ((fabs(PrimDirnCosn.YUnit.X - xnorm) <= 1.0e-3) &&
859 (fabs(PrimDirnCosn.YUnit.Y - ynorm) <= 1.0e-3) &&
860 (fabs(PrimDirnCosn.YUnit.Z - znorm) <= 1.0e-3))
861 flagDC = 1;
862 if (DebugLevel == 202) {
863 printf("First attempt: \n");
864 PrintDirnCosn3D(PrimDirnCosn);
865 printf("\n");
866 }
867
868 if (!flagDC) // if DC search is unsuccessful, try the other orientation
869 {
870 SurfLX = sqrt((xvert[2] - xvert[1]) * (xvert[2] - xvert[1]) +
871 (yvert[2] - yvert[1]) * (yvert[2] - yvert[1]) +
872 (zvert[2] - zvert[1]) * (zvert[2] - zvert[1]));
873 SurfLZ = sqrt((xvert[0] - xvert[1]) * (xvert[0] - xvert[1]) +
874 (yvert[0] - yvert[1]) * (yvert[0] - yvert[1]) +
875 (zvert[0] - zvert[1]) * (zvert[0] - zvert[1]));
876 // Direction cosines - note that right angle corner is [1]-th element
877 PrimDirnCosn.XUnit.X = (xvert[2] - xvert[1]) / SurfLX;
878 PrimDirnCosn.XUnit.Y = (yvert[2] - yvert[1]) / SurfLX;
879 PrimDirnCosn.XUnit.Z = (zvert[2] - zvert[1]) / SurfLX;
880 PrimDirnCosn.ZUnit.X = (xvert[0] - xvert[1]) / SurfLZ;
881 PrimDirnCosn.ZUnit.Y = (yvert[0] - yvert[1]) / SurfLZ;
882 PrimDirnCosn.ZUnit.Z = (zvert[0] - zvert[1]) / SurfLZ;
883 PrimDirnCosn.YUnit =
884 Vector3DCrossProduct(&PrimDirnCosn.ZUnit, &PrimDirnCosn.XUnit);
885 if ((fabs(PrimDirnCosn.YUnit.X - xnorm) <= 1.0e-3) &&
886 (fabs(PrimDirnCosn.YUnit.Y - ynorm) <= 1.0e-3) &&
887 (fabs(PrimDirnCosn.YUnit.Z - znorm) <= 1.0e-3))
888 flagDC = 2;
889 if (DebugLevel == 202) {
890 printf("Second attempt: \n");
891 PrintDirnCosn3D(PrimDirnCosn);
892 printf("\n");
893 }
894 }
895
896 if (!flagDC) // No other possibility, DC search failed!!!
897 {
898 printf("Triangle DC problem ... returning ...\n");
899 // getchar();
900 return -1;
901 }
902
903 // primitive direction cosine assignments
904 PrimDC[prim].XUnit.X = PrimDirnCosn.XUnit.X;
905 PrimDC[prim].XUnit.Y = PrimDirnCosn.XUnit.Y;
906 PrimDC[prim].XUnit.Z = PrimDirnCosn.XUnit.Z;
907 PrimDC[prim].YUnit.X = PrimDirnCosn.YUnit.X;
908 PrimDC[prim].YUnit.Y = PrimDirnCosn.YUnit.Y;
909 PrimDC[prim].YUnit.Z = PrimDirnCosn.YUnit.Z;
910 PrimDC[prim].ZUnit.X = PrimDirnCosn.ZUnit.X;
911 PrimDC[prim].ZUnit.Y = PrimDirnCosn.ZUnit.Y;
912 PrimDC[prim].ZUnit.Z = PrimDirnCosn.ZUnit.Z;
913
914 // primitive origin - for a triangle, origin is at the right corner
915 PrimOriginX[prim] = SurfX;
916 PrimOriginY[prim] = SurfY;
917 PrimOriginZ[prim] = SurfZ;
918 PrimLX[prim] = SurfLX;
919 PrimLZ[prim] = SurfLZ;
920 if (flagDC == 1) {
921 int tmpVolRef1 = VolRef1[prim];
922 VolRef1[prim] = VolRef2[prim];
923 VolRef2[prim] = tmpVolRef1;
924 int tmpEpsilon1 = Epsilon1[prim];
925 Epsilon1[prim] = Epsilon2[prim];
926 Epsilon2[prim] = tmpEpsilon1;
927 }
928
929 double SurfLambda = lambda;
930 double SurfV = potential;
931
932 // file output for a primitive
933 if (OptPrimitiveFiles) {
934 char OutPrim[256];
935 strcpy(OutPrim, ModelOutDir);
936 strcat(OutPrim, "/Primitives/Primitive");
937 strcat(OutPrim, primstr);
938 strcat(OutPrim, ".out");
939 fPrim = fopen(OutPrim, "w");
940 if (fPrim == NULL) {
941 neBEMMessage("DiscretizeTriangle - OutPrim");
942 return -1;
943 }
944 fprintf(fPrim, "#prim: %d, nvertex: %d\n", prim, nvertex);
945 fprintf(fPrim, "Node1: %lg\t%lg\t%lg\n", xvert[0], yvert[0], zvert[0]);
946 fprintf(fPrim, "Node2: %lg\t%lg\t%lg\n", xvert[1], yvert[1], zvert[1]);
947 fprintf(fPrim, "Node3: %lg\t%lg\t%lg\n", xvert[2], yvert[2], zvert[2]);
948 fprintf(fPrim, "PrimOrigin: %lg\t%lg\t%lg\n", PrimOriginX[prim],
949 PrimOriginY[prim], PrimOriginZ[prim]);
950 fprintf(fPrim, "Primitive lengths: %lg\t%lg\n", PrimLX[prim], PrimLZ[prim]);
951 fprintf(fPrim, "Norm: %lg\t%lg\t%lg\n", xnorm, ynorm, znorm);
952 fprintf(fPrim, "#volref1: %d, volref2: %d\n", volref1, volref2);
953 fprintf(fPrim, "#NbSegX: %d, NbSegZ: %d (check note!)\n", NbSegX, NbSegZ);
954 fprintf(fPrim, "#ParentObj: %d\tEType: %d\n", SurfParentObj, SurfEType);
955 fprintf(fPrim, "#SurfX\tSurfY\tSurfZ\tSurfLX\tSurfLZ (Rt. Corner)\n");
956 fprintf(fPrim, "%lg\t%lg\t%lg\t%lg\t%lg\n", SurfX, SurfY, SurfZ, SurfLX,
957 SurfLZ);
958 fprintf(fPrim, "#DirnCosn: \n");
959 fprintf(fPrim, "%lg, %lg, %lg\n", PrimDirnCosn.XUnit.X,
960 PrimDirnCosn.XUnit.Y, PrimDirnCosn.XUnit.Z);
961 fprintf(fPrim, "%lg, %lg, %lg\n", PrimDirnCosn.YUnit.X,
962 PrimDirnCosn.YUnit.Y, PrimDirnCosn.YUnit.Z);
963 fprintf(fPrim, "%lg, %lg, %lg\n", PrimDirnCosn.ZUnit.X,
964 PrimDirnCosn.ZUnit.Y, PrimDirnCosn.ZUnit.Z);
965 fprintf(fPrim, "#SurfLambda: %lg\tSurfV: %lg\n", SurfLambda, SurfV);
966 }
967
968 // necessary for gnuplot
969 if (OptGnuplot && OptGnuplotPrimitives) {
970 char gpPrim[256];
971 strcpy(gpPrim, MeshOutDir);
972 strcat(gpPrim, "/GViewDir/gpPrim");
973 strcat(gpPrim, primstr);
974 strcat(gpPrim, ".out");
975 fgpPrim = fopen(gpPrim, "w");
976 if (fgpPrim == NULL) {
977 neBEMMessage("DiscretizeTriangle - OutgpPrim");
978 return -1;
979 }
980 fprintf(fgpPrim, "%g\t%g\t%g\n\n", xvert[0], yvert[0], zvert[0]);
981 fprintf(fgpPrim, "%g\t%g\t%g\n\n", xvert[1], yvert[1], zvert[1]);
982 fprintf(fgpPrim, "%g\t%g\t%g\n\n", xvert[2], yvert[2], zvert[2]);
983 fprintf(fgpPrim, "%g\t%g\t%g\n", xvert[0], yvert[0], zvert[0]);
984 fclose(fgpPrim);
985
986 if (prim == 1)
987 fprintf(fgnuPrim, " '%s\' w l", gpPrim);
988 else
989 fprintf(fgnuPrim, ", \\\n \'%s\' w l", gpPrim);
990 }
991
992 // file outputs for elements on primitive
993 if (OptElementFiles) {
994 char OutElem[256];
995 strcpy(OutElem, MeshOutDir);
996 strcat(OutElem, "/Elements/ElemOnPrim");
997 strcat(OutElem, primstr);
998 strcat(OutElem, ".out");
999 fElem = fopen(OutElem, "w");
1000 if (fElem == NULL) {
1001 neBEMMessage("DiscretizeTriangle - OutElem");
1002 return -1;
1003 }
1004 }
1005 // gnuplot friendly file outputs for elements on primitive
1006 if (OptGnuplot && OptGnuplotElements) {
1007 strcpy(gpElem, MeshOutDir);
1008 strcat(gpElem, "/GViewDir/gpElemOnPrim");
1009 strcat(gpElem, primstr);
1010 strcat(gpElem, ".out");
1011 fgpElem = fopen(gpElem, "w");
1012 // assert(fgpElem != NULL);
1013 if (fgpElem == NULL) {
1014 neBEMMessage("DiscretizeTriangle - OutgpElem");
1015 if (fElem) fclose(fElem);
1016 return -1;
1017 }
1018 // gnuplot friendly file outputs for elements on primitive
1019 strcpy(gpMesh, MeshOutDir);
1020 strcat(gpMesh, "/GViewDir/gpMeshOnPrim");
1021 strcat(gpMesh, primstr);
1022 strcat(gpMesh, ".out");
1023 fgpMesh = fopen(gpMesh, "w");
1024 if (fgpMesh == NULL) {
1025 neBEMMessage("DiscretizeTriangle - OutgpMesh");
1026 fclose(fgpElem);
1027 if (fElem) fclose(fElem);
1028 return -1;
1029 }
1030 }
1031
1032 // Compute element positions (CGs) in primitive local coordinate system (PCS).
1033 // Then map these CGs to the global coordinate system.
1034 // (xav, 0, zav) is the CG of an element wrt the primitive coordinate system
1035 // From this, we find the offset of the element from the primitive CG in the
1036 // global coordinate system by just rotating the displacement vector
1037 // (0,0,0) -> (xav,0,zav) using the known DCM of the surface.
1038 // This rotated vector (now in the global system) is then added to the
1039 // position vector of the primitive CG (always in the global system) to get
1040 // the CG of the element in the global system.
1041
1042 // Consult note for clarifications on the discretization algorithm.
1043 // SurfElLX is true for the lowest row of elements, but indicative of the
1044 // other rows
1045 // Make sure that the larger number is being used to discretize the longer
1046 // side - can be OverSmart in certain cases - there may be cases where
1047 // smaller number of elements suffice for a longer side
1048 if (NbSegX == NbSegZ) {
1049 SurfElLX = SurfLX / NbSegX; // element sizes
1050 SurfElLZ = SurfLZ / NbSegZ;
1051 } else if (NbSegX > NbSegZ) {
1052 if (SurfLX > SurfLZ) {
1053 SurfElLX = SurfLX / NbSegX; // element sizes
1054 SurfElLZ = SurfLZ / NbSegZ;
1055 } else // interchange NbSegX and NbSegZ
1056 {
1057 int tmp = NbSegZ;
1058 NbSegZ = NbSegX;
1059 NbSegX = tmp;
1060 SurfElLX = SurfLX / NbSegX; // element sizes
1061 SurfElLZ = SurfLZ / NbSegZ;
1062 }
1063 } // NbSegX > NbSegZ
1064 else // NbSegX < NbSegZ
1065 {
1066 if (SurfLX < SurfLZ) {
1067 SurfElLX = SurfLX / NbSegX; // element sizes
1068 SurfElLZ = SurfLZ / NbSegZ;
1069 } else // interchange NbSegX and NbSegZ
1070 {
1071 int tmp = NbSegZ;
1072 NbSegZ = NbSegX;
1073 NbSegX = tmp;
1074 SurfElLX = SurfLX / NbSegX; // element sizes
1075 SurfElLZ = SurfLZ / NbSegZ;
1076 }
1077 }
1078
1079 // Analysis of element aspect ratio.
1080 // Note that we can afford only to reduce the total number of elements.
1081 // Otherwise, we'll have to realloc `EleArr' array.
1082 // Later, we'll make the corrections such that the total number of elements
1083 // remain close to the originally intended value.
1084 double AR = SurfElLX / SurfElLZ; // indicative element aspect ratio
1085 if (OptPrimitiveFiles) {
1086 fprintf(fPrim,
1087 "Using the input, the aspect ratio of the elements on prim: %d\n",
1088 prim);
1089 fprintf(fPrim,
1090 "NbSegX: %d, SurfElLX: %lg, NbSegZ: %d, SurfElLZ: %lg, AR: %lg\n",
1091 NbSegX, SurfElLX, NbSegZ, SurfElLZ, AR);
1092 }
1093 if (AR > 10.0) // NbSegZ needs to be reduced
1094 {
1095 double tmpElLZ = SurfElLX / 10.0;
1096 NbSegZ = (int)(SurfLZ / tmpElLZ);
1097 if (NbSegZ <= 0) NbSegZ = 1;
1098 SurfElLZ = SurfLZ / NbSegZ;
1099 }
1100 if (AR < 0.1) // NbSegX need to be reduced
1101 {
1102 double tmpElLX = SurfElLZ * 0.1;
1103 NbSegX = (int)(SurfLX / tmpElLX);
1104 if (NbSegX <= 0) NbSegX = 1;
1105 SurfElLX = SurfLX / NbSegX;
1106 }
1107 AR = SurfElLX / SurfElLZ; // indicative element aspect ratio
1108 if (OptPrimitiveFiles) {
1109 fprintf(fPrim, "After analyzing the likely aspect ratio of the elements\n");
1110 fprintf(fPrim,
1111 "NbSegX: %d, SurfElLX: %lg, NbSegZ: %d, SurfElLZ: %lg, AR: %lg\n",
1112 NbSegX, SurfElLX, NbSegZ, SurfElLZ, AR);
1113 }
1114
1115 // The top-most right angle triangle is a lone element on the highest row, its
1116 // properties being determined irrespective of the others. Despite that, this
1117 // element is being treated as one among the others, especially to facilitate
1118 // element indexing, file output etc.
1119 // Each row, thus, has at least one triangle, and possibly several rectangular
1120 // elements. All the elements in any given row has the same SurfElLZ as has
1121 // been determined above.
1122 // First we create the triangular element and then move on to create the
1123 // rectangular ones.
1124 double xv0, yv0, zv0, xv1, yv1, zv1, xv2, yv2, zv2;
1125 ElementBgn[prim] = EleCntr + 1;
1126 for (int k = 1; k <= NbSegZ; ++k) // consider the k-th row
1127 {
1128 double grad = (SurfLZ / SurfLX);
1129 double zlopt = (k - 1) * SurfElLZ;
1130 double zhipt = (k)*SurfElLZ;
1131 double xlopt = (SurfLZ - zlopt) / grad; // used again on 21 Feb 2014
1132 double xhipt = (SurfLZ - zhipt) / grad;
1133
1134 // the triangular element on the k-th row can now be specified in PCS
1135 double xtorigin = xhipt;
1136 double ytorigin = 0.0;
1137 double ztorigin = zlopt;
1138 { // Separate block for position rotation - local2global
1139 Point3D localDisp, globalDisp;
1140
1141 localDisp.X = xtorigin;
1142 localDisp.Y = ytorigin;
1143 localDisp.Z = ztorigin;
1144 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1145 SurfElX =
1146 SurfX + globalDisp.X; // these are the coords in GCS of origin of
1147 SurfElY =
1148 SurfY + globalDisp.Y; // the triangluar element under consideration
1149 SurfElZ = SurfZ + globalDisp.Z;
1150 } // vector rotation over
1151
1152 // Assign element values and write in the file
1153 ++EleCntr;
1154 if (EleCntr > NbElements) {
1155 neBEMMessage("DiscretizeTriangle - EleCntr more than NbElements 1!");
1156 if (fgpElem) fclose(fgpElem);
1157 if (fgpMesh) fclose(fgpMesh);
1158 return -1;
1159 }
1160
1161 (EleArr + EleCntr - 1)->DeviceNb =
1162 1; // At present, there is only one device
1163 (EleArr + EleCntr - 1)->ComponentNb = SurfParentObj;
1164 (EleArr + EleCntr - 1)->PrimitiveNb = prim;
1165 (EleArr + EleCntr - 1)->Id = EleCntr;
1166 (EleArr + EleCntr - 1)->G.Type = 3; // triangular here
1167 (EleArr + EleCntr - 1)->G.Origin.X = SurfElX;
1168 (EleArr + EleCntr - 1)->G.Origin.Y = SurfElY;
1169 (EleArr + EleCntr - 1)->G.Origin.Z = SurfElZ;
1170 // (EleArr+EleCntr-1)->G.LX = SurfElLX; // previously written as xlopt -
1171 // xhipt;
1172 (EleArr + EleCntr - 1)->G.LX =
1173 xlopt - xhipt; // back to old ways on 21 Feb 2014
1174 (EleArr + EleCntr - 1)->G.LZ = SurfElLZ;
1175 (EleArr + EleCntr - 1)->G.LZ =
1176 zhipt - zlopt; // to be on the safe side, 21/2/14
1177 (EleArr + EleCntr - 1)->G.dA =
1178 0.5 * (EleArr + EleCntr - 1)->G.LX * (EleArr + EleCntr - 1)->G.LZ;
1179 // Safe to use the direction cosines obtained for the triangular primitive
1180 // since they are bound to remain unchanged for the rectangular sub-elements
1181 (EleArr + EleCntr - 1)->G.DC.XUnit.X = PrimDirnCosn.XUnit.X;
1182 (EleArr + EleCntr - 1)->G.DC.XUnit.Y = PrimDirnCosn.XUnit.Y;
1183 (EleArr + EleCntr - 1)->G.DC.XUnit.Z = PrimDirnCosn.XUnit.Z;
1184 (EleArr + EleCntr - 1)->G.DC.YUnit.X = PrimDirnCosn.YUnit.X;
1185 (EleArr + EleCntr - 1)->G.DC.YUnit.Y = PrimDirnCosn.YUnit.Y;
1186 (EleArr + EleCntr - 1)->G.DC.YUnit.Z = PrimDirnCosn.YUnit.Z;
1187 (EleArr + EleCntr - 1)->G.DC.ZUnit.X = PrimDirnCosn.ZUnit.X;
1188 (EleArr + EleCntr - 1)->G.DC.ZUnit.Y = PrimDirnCosn.ZUnit.Y;
1189 (EleArr + EleCntr - 1)->G.DC.ZUnit.Z = PrimDirnCosn.ZUnit.Z;
1190 (EleArr + EleCntr - 1)->E.Type = SurfEType;
1191 (EleArr + EleCntr - 1)->E.Lambda = SurfLambda;
1192 (EleArr + EleCntr - 1)->Solution = 0.0;
1193 (EleArr + EleCntr - 1)->Assigned = charge;
1194 // Boundary condition is applied at the barycenter, not at the origin
1195 // of the element coordinate system (ECS) which is at the right corner
1196 (EleArr + EleCntr - 1)->BC.NbOfBCs = 1; // assume one BC per element
1197
1198 xv0 = (EleArr + EleCntr - 1)->G.Origin.X;
1199 yv0 = (EleArr + EleCntr - 1)->G.Origin.Y;
1200 zv0 = (EleArr + EleCntr - 1)->G.Origin.Z;
1201 xv1 = (EleArr + EleCntr - 1)->G.Origin.X +
1202 (EleArr + EleCntr - 1)->G.LX * (EleArr + EleCntr - 1)->G.DC.XUnit.X;
1203 yv1 = (EleArr + EleCntr - 1)->G.Origin.Y +
1204 (EleArr + EleCntr - 1)->G.LX * (EleArr + EleCntr - 1)->G.DC.XUnit.Y;
1205 zv1 = (EleArr + EleCntr - 1)->G.Origin.Z +
1206 (EleArr + EleCntr - 1)->G.LX * (EleArr + EleCntr - 1)->G.DC.XUnit.Z;
1207 xv2 = (EleArr + EleCntr - 1)->G.Origin.X +
1208 (EleArr + EleCntr - 1)->G.LZ * (EleArr + EleCntr - 1)->G.DC.ZUnit.X;
1209 yv2 = (EleArr + EleCntr - 1)->G.Origin.Y +
1210 (EleArr + EleCntr - 1)->G.LZ * (EleArr + EleCntr - 1)->G.DC.ZUnit.Y;
1211 zv2 = (EleArr + EleCntr - 1)->G.Origin.Z +
1212 (EleArr + EleCntr - 1)->G.LZ * (EleArr + EleCntr - 1)->G.DC.ZUnit.Z;
1213 // assign vertices of the element
1214 (EleArr + EleCntr - 1)->G.Vertex[0].X = xv0;
1215 (EleArr + EleCntr - 1)->G.Vertex[0].Y = yv0;
1216 (EleArr + EleCntr - 1)->G.Vertex[0].Z = zv0;
1217 (EleArr + EleCntr - 1)->G.Vertex[1].X = xv1;
1218 (EleArr + EleCntr - 1)->G.Vertex[1].Y = yv1;
1219 (EleArr + EleCntr - 1)->G.Vertex[1].Z = zv1;
1220 (EleArr + EleCntr - 1)->G.Vertex[2].X = xv2;
1221 (EleArr + EleCntr - 1)->G.Vertex[2].Y = yv2;
1222 (EleArr + EleCntr - 1)->G.Vertex[2].Z = zv2;
1223 (EleArr + EleCntr - 1)->G.Vertex[3].X = 0.0;
1224 (EleArr + EleCntr - 1)->G.Vertex[3].Y = 0.0;
1225 (EleArr + EleCntr - 1)->G.Vertex[3].Z = 0.0;
1226
1227 if (DebugLevel == 201) {
1228 printf("Primitive nb: %d\n", (EleArr + EleCntr - 1)->PrimitiveNb);
1229 printf("Element id: %d\n", (EleArr + EleCntr - 1)->Id);
1230 printf("Element X, Y, Z: %lg %lg %lg\n",
1231 (EleArr + EleCntr - 1)->G.Origin.X,
1232 (EleArr + EleCntr - 1)->G.Origin.Y,
1233 (EleArr + EleCntr - 1)->G.Origin.Z);
1234 printf("Element LX, LZ: %lg %lg\n", (EleArr + EleCntr - 1)->G.LX,
1235 (EleArr + EleCntr - 1)->G.LZ);
1236 printf("Element (primitive) X axis dirn cosines: %lg, %lg, %lg\n",
1237 PrimDirnCosn.XUnit.X, PrimDirnCosn.XUnit.Y, PrimDirnCosn.XUnit.Z);
1238 printf("Element (primitive) Y axis dirn cosines: %lg, %lg, %lg\n",
1239 PrimDirnCosn.YUnit.X, PrimDirnCosn.YUnit.Y, PrimDirnCosn.YUnit.Z);
1240 printf("Element (primitive) Z axis dirn cosines: %lg, %lg, %lg\n",
1241 PrimDirnCosn.ZUnit.X, PrimDirnCosn.ZUnit.Y, PrimDirnCosn.ZUnit.Z);
1242 }
1243 // Following are the location in the ECS
1244 double dxl = (EleArr + EleCntr - 1)->G.LX / 3.0;
1245 double dyl = 0.0;
1246 double dzl = (EleArr + EleCntr - 1)->G.LZ / 3.0;
1247 { // Separate block for position rotation - local2global
1248 Point3D localDisp, globalDisp;
1249
1250 localDisp.X = dxl;
1251 localDisp.Y = dyl;
1252 localDisp.Z = dzl;
1253 if (DebugLevel == 201) {
1254 printf("Element dxl, dxy, dxz: %lg %lg %lg\n", localDisp.X, localDisp.Y,
1255 localDisp.Z);
1256 }
1257
1258 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1259 (EleArr + EleCntr - 1)->BC.CollPt.X =
1260 (EleArr + EleCntr - 1)->G.Origin.X + globalDisp.X;
1261 (EleArr + EleCntr - 1)->BC.CollPt.Y =
1262 (EleArr + EleCntr - 1)->G.Origin.Y + globalDisp.Y;
1263 (EleArr + EleCntr - 1)->BC.CollPt.Z =
1264 (EleArr + EleCntr - 1)->G.Origin.Z + globalDisp.Z;
1265 if (DebugLevel == 201) {
1266 printf("Element global dxl, dxy, dxz: %lg %lg %lg\n", globalDisp.X,
1267 globalDisp.Y, globalDisp.Z);
1268 printf("Element BCX, BCY, BCZ: %lg %lg %lg\n",
1269 (EleArr + EleCntr - 1)->BC.CollPt.X,
1270 (EleArr + EleCntr - 1)->BC.CollPt.Y,
1271 (EleArr + EleCntr - 1)->BC.CollPt.Z);
1272 }
1273 } // vector rotation over
1274 // (EleArr+EleCntr-1)->BC.Value = SurfV; // assigned in BoundaryConditions
1275
1276 if (OptElementFiles) {
1277 fprintf(fElem, "##Element Counter: %d\n", EleCntr);
1278 fprintf(fElem, "#DevNb\tCompNb\tPrimNb\tId\n");
1279 fprintf(fElem, "%d\t%d\t%d\t%d\n", (EleArr + EleCntr - 1)->DeviceNb,
1280 (EleArr + EleCntr - 1)->ComponentNb,
1281 (EleArr + EleCntr - 1)->PrimitiveNb, (EleArr + EleCntr - 1)->Id);
1282 fprintf(fElem, "#GType\tX\tY\tZ\tLX\tLZ\tdA\n");
1283 fprintf(fElem, "%d\t%.16lg\t%.16lg\t%.16lg\t%.16lg\t%.16lg\t%.16lg\n",
1284 (EleArr + EleCntr - 1)->G.Type,
1285 (EleArr + EleCntr - 1)->G.Origin.X,
1286 (EleArr + EleCntr - 1)->G.Origin.Y,
1287 (EleArr + EleCntr - 1)->G.Origin.Z, (EleArr + EleCntr - 1)->G.LX,
1288 (EleArr + EleCntr - 1)->G.LZ, (EleArr + EleCntr - 1)->G.dA);
1289 fprintf(fElem, "#DirnCosn: \n");
1290 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.XUnit.X,
1291 (EleArr + EleCntr - 1)->G.DC.XUnit.Y,
1292 (EleArr + EleCntr - 1)->G.DC.XUnit.Z);
1293 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.YUnit.X,
1294 (EleArr + EleCntr - 1)->G.DC.YUnit.Y,
1295 (EleArr + EleCntr - 1)->G.DC.YUnit.Z);
1296 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.ZUnit.X,
1297 (EleArr + EleCntr - 1)->G.DC.ZUnit.Y,
1298 (EleArr + EleCntr - 1)->G.DC.ZUnit.Z);
1299 fprintf(fElem, "#EType\tLambda\n");
1300 fprintf(fElem, "%d\t%lg\n", (EleArr + EleCntr - 1)->E.Type,
1301 (EleArr + EleCntr - 1)->E.Lambda);
1302 fprintf(fElem, "#NbBCs\tCPX\tCPY\tCPZ\tValue\n");
1303 fprintf(fElem, "%d\t%.16lg\t%.16lg\t%.16lg\t%lg\n",
1304 (EleArr + EleCntr - 1)->BC.NbOfBCs,
1305 (EleArr + EleCntr - 1)->BC.CollPt.X,
1306 (EleArr + EleCntr - 1)->BC.CollPt.Y,
1307 (EleArr + EleCntr - 1)->BC.CollPt.Z,
1308 (EleArr + EleCntr - 1)->BC.Value);
1309 } // if OptElementFiles
1310
1311 // mark bary-center and draw mesh
1312 if (OptGnuplot && OptGnuplotElements) {
1313 fprintf(fgpElem, "%g\t%g\t%g\n", (EleArr + EleCntr - 1)->BC.CollPt.X,
1314 (EleArr + EleCntr - 1)->BC.CollPt.Y,
1315 (EleArr + EleCntr - 1)->BC.CollPt.Z);
1316
1317 // draw mesh
1318 // assign vertices of the element
1319 xv0 = (EleArr + EleCntr - 1)->G.Vertex[0].X;
1320 yv0 = (EleArr + EleCntr - 1)->G.Vertex[0].Y;
1321 zv0 = (EleArr + EleCntr - 1)->G.Vertex[0].Z;
1322 xv1 = (EleArr + EleCntr - 1)->G.Vertex[1].X;
1323 yv1 = (EleArr + EleCntr - 1)->G.Vertex[1].Y;
1324 zv1 = (EleArr + EleCntr - 1)->G.Vertex[1].Z;
1325 xv2 = (EleArr + EleCntr - 1)->G.Vertex[2].X;
1326 yv2 = (EleArr + EleCntr - 1)->G.Vertex[2].Y;
1327 zv2 = (EleArr + EleCntr - 1)->G.Vertex[2].Z;
1328
1329 fprintf(fgpMesh, "%g\t%g\t%g\n", xv0, yv0, zv0);
1330 fprintf(fgpMesh, "%g\t%g\t%g\n", xv1, yv1, zv1);
1331 fprintf(fgpMesh, "%g\t%g\t%g\n", xv2, yv2, zv2);
1332 fprintf(fgpMesh, "%g\t%g\t%g\n\n", xv0, yv0, zv0);
1333 } // if OptGnuplotElements
1334
1335 if (k == NbSegZ) // no rectangular element on this row
1336 continue;
1337
1338 // determine NbSegXOnThisRow and ElLXOnThisRow for the rectagnular elements
1339 // and then loop.
1340 double RowLX = xhipt; // the triangular portion is left outside the slice
1341 int NbSegXOnThisRow;
1342 double ElLXOnThisRow;
1343 if (RowLX <= SurfElLX) {
1344 NbSegXOnThisRow = 1;
1345 ElLXOnThisRow = RowLX;
1346 } else {
1347 NbSegXOnThisRow = (int)(RowLX / SurfElLX);
1348 ElLXOnThisRow = RowLX / NbSegXOnThisRow;
1349 }
1350 double x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3;
1351 for (int i = 1; i <= NbSegXOnThisRow; ++i) {
1352 double xorigin = (i - 1) * ElLXOnThisRow + 0.5 * ElLXOnThisRow; // PCS
1353 double yorigin = 0.0; // centroid of the rectagnular element
1354 double zorigin = 0.5 * (zlopt + zhipt);
1355 // printf("k: %d, i: %d, xo: %lg, yo: %lg, zo: %lg\n", k, i,
1356 // xorigin, yorigin, zorigin);
1357 // printf("xlopt: %lg, zlopt: %lg, xhipt: %lg, zhipt: %lg\n",
1358 // xlopt, zlopt, xhipt, zhipt);
1359
1360 { // Separate block for vector rotation
1361 Point3D localDisp, globalDisp;
1362
1363 localDisp.X = xorigin;
1364 localDisp.Y = yorigin;
1365 localDisp.Z = zorigin;
1366 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1367 SurfElX = SurfX + globalDisp.X; // GCS
1368 SurfElY = SurfY + globalDisp.Y;
1369 SurfElZ = SurfZ + globalDisp.Z;
1370 } // vector rotation over
1371 // printf("SurfX: %lg, SurfY: %lg, SurfZ: %lg\n",
1372 // SurfX, SurfY, SurfZ);
1373 // printf("SurfElX: %lg, SurfElY: %lg, SurfElZ: %lg\n",
1374 // SurfElX, SurfElY, SurfElZ);
1375
1376 // Assign element values and write in the file
1377 ++EleCntr;
1378 if (EleCntr > NbElements) {
1379 neBEMMessage("DiscretizeTriangle - EleCntr more than NbElements 2!");
1380 return -1;
1381 }
1382
1383 (EleArr + EleCntr - 1)->DeviceNb =
1384 1; // At present, there is only one device
1385 (EleArr + EleCntr - 1)->ComponentNb = SurfParentObj;
1386 (EleArr + EleCntr - 1)->PrimitiveNb = prim;
1387 (EleArr + EleCntr - 1)->Id = EleCntr;
1388 (EleArr + EleCntr - 1)->G.Type = 4; // rectagnular here
1389 (EleArr + EleCntr - 1)->G.Origin.X = SurfElX;
1390 (EleArr + EleCntr - 1)->G.Origin.Y = SurfElY;
1391 (EleArr + EleCntr - 1)->G.Origin.Z = SurfElZ;
1392 (EleArr + EleCntr - 1)->G.LX = ElLXOnThisRow;
1393 (EleArr + EleCntr - 1)->G.LZ = SurfElLZ;
1394 (EleArr + EleCntr - 1)->G.LZ =
1395 zhipt - zlopt; // to be on the safe side! 21/2/14
1396 (EleArr + EleCntr - 1)->G.dA =
1397 (EleArr + EleCntr - 1)->G.LX * (EleArr + EleCntr - 1)->G.LZ;
1398 // Safe to use the direction cosines obtained for the triangular primitive
1399 // since they are bound to remain unchanged for the triangular
1400 // sub-elements
1401 (EleArr + EleCntr - 1)->G.DC.XUnit.X = PrimDirnCosn.XUnit.X;
1402 (EleArr + EleCntr - 1)->G.DC.XUnit.Y = PrimDirnCosn.XUnit.Y;
1403 (EleArr + EleCntr - 1)->G.DC.XUnit.Z = PrimDirnCosn.XUnit.Z;
1404 (EleArr + EleCntr - 1)->G.DC.YUnit.X = PrimDirnCosn.YUnit.X;
1405 (EleArr + EleCntr - 1)->G.DC.YUnit.Y = PrimDirnCosn.YUnit.Y;
1406 (EleArr + EleCntr - 1)->G.DC.YUnit.Z = PrimDirnCosn.YUnit.Z;
1407 (EleArr + EleCntr - 1)->G.DC.ZUnit.X = PrimDirnCosn.ZUnit.X;
1408 (EleArr + EleCntr - 1)->G.DC.ZUnit.Y = PrimDirnCosn.ZUnit.Y;
1409 (EleArr + EleCntr - 1)->G.DC.ZUnit.Z = PrimDirnCosn.ZUnit.Z;
1410 (EleArr + EleCntr - 1)->E.Type = SurfEType;
1411 (EleArr + EleCntr - 1)->E.Lambda = SurfLambda;
1412 (EleArr + EleCntr - 1)->Solution = 0.0;
1413 (EleArr + EleCntr - 1)->Assigned = charge;
1414 // Boundary condition is applied at the origin for this rectangular
1415 // element coordinate system (ECS)
1416 (EleArr + EleCntr - 1)->BC.NbOfBCs = 1; // assume one BC per element
1417 // Following are the location in the ECS
1418 (EleArr + EleCntr - 1)->BC.CollPt.X = (EleArr + EleCntr - 1)->G.Origin.X;
1419 (EleArr + EleCntr - 1)->BC.CollPt.Y = (EleArr + EleCntr - 1)->G.Origin.Y;
1420 (EleArr + EleCntr - 1)->BC.CollPt.Z = (EleArr + EleCntr - 1)->G.Origin.Z;
1421 // find element vertices
1422 // 1) displacement vector in the ECS is first identified
1423 // 2) this vector is transformed to the GCS
1424 // 3) the global displacement vector, when added to the centroid in GCS,
1425 // gives the node positions in GCS.
1426 x0 = -0.5 *
1427 (EleArr + EleCntr - 1)->G.LX; // xyz displacement of node wrt ECS
1428 y0 = 0.0;
1429 z0 = -0.5 * (EleArr + EleCntr - 1)->G.LZ;
1430 { // Separate block for position rotation - local2global
1431 Point3D localDisp, globalDisp;
1432
1433 localDisp.X = x0;
1434 localDisp.Y = y0;
1435 localDisp.Z = z0; // displacement in GCS
1436 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1437 x0 = (EleArr + EleCntr - 1)->G.Origin.X +
1438 globalDisp.X; // xyz position in GCS
1439 y0 = (EleArr + EleCntr - 1)->G.Origin.Y + globalDisp.Y;
1440 z0 = (EleArr + EleCntr - 1)->G.Origin.Z + globalDisp.Z;
1441 } // position rotation over
1442
1443 x1 = 0.5 * (EleArr + EleCntr - 1)->G.LX;
1444 y1 = 0.0;
1445 z1 = -0.5 * (EleArr + EleCntr - 1)->G.LZ;
1446 { // Separate block for position rotation - local2global
1447 Point3D localDisp, globalDisp;
1448
1449 localDisp.X = x1;
1450 localDisp.Y = y1;
1451 localDisp.Z = z1;
1452 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1453 x1 = (EleArr + EleCntr - 1)->G.Origin.X + globalDisp.X;
1454 y1 = (EleArr + EleCntr - 1)->G.Origin.Y + globalDisp.Y;
1455 z1 = (EleArr + EleCntr - 1)->G.Origin.Z + globalDisp.Z;
1456 } // position rotation over
1457
1458 x2 = 0.5 * (EleArr + EleCntr - 1)->G.LX;
1459 y2 = 0.0;
1460 z2 = 0.5 * (EleArr + EleCntr - 1)->G.LZ;
1461 { // Separate block for position rotation - local2global
1462 Point3D localDisp, globalDisp;
1463
1464 localDisp.X = x2;
1465 localDisp.Y = y2;
1466 localDisp.Z = z2;
1467 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1468 x2 = (EleArr + EleCntr - 1)->G.Origin.X + globalDisp.X;
1469 y2 = (EleArr + EleCntr - 1)->G.Origin.Y + globalDisp.Y;
1470 z2 = (EleArr + EleCntr - 1)->G.Origin.Z + globalDisp.Z;
1471 } // position rotation over
1472
1473 x3 = -0.5 * (EleArr + EleCntr - 1)->G.LX;
1474 y3 = 0.0;
1475 z3 = 0.5 * (EleArr + EleCntr - 1)->G.LZ;
1476 { // Separate block for position rotation - local2global
1477 Point3D localDisp, globalDisp;
1478
1479 localDisp.X = x3;
1480 localDisp.Y = y3;
1481 localDisp.Z = z3;
1482 globalDisp = RotatePoint3D(&localDisp, &PrimDirnCosn, local2global);
1483 x3 = (EleArr + EleCntr - 1)->G.Origin.X + globalDisp.X;
1484 y3 = (EleArr + EleCntr - 1)->G.Origin.Y + globalDisp.Y;
1485 z3 = (EleArr + EleCntr - 1)->G.Origin.Z + globalDisp.Z;
1486 } // position rotation over
1487
1488 // assign vertices of the element
1489 (EleArr + EleCntr - 1)->G.Vertex[0].X = x0;
1490 (EleArr + EleCntr - 1)->G.Vertex[0].Y = y0;
1491 (EleArr + EleCntr - 1)->G.Vertex[0].Z = z0;
1492 (EleArr + EleCntr - 1)->G.Vertex[1].X = x1;
1493 (EleArr + EleCntr - 1)->G.Vertex[1].Y = y1;
1494 (EleArr + EleCntr - 1)->G.Vertex[1].Z = z1;
1495 (EleArr + EleCntr - 1)->G.Vertex[2].X = x2;
1496 (EleArr + EleCntr - 1)->G.Vertex[2].Y = y2;
1497 (EleArr + EleCntr - 1)->G.Vertex[2].Z = z2;
1498 (EleArr + EleCntr - 1)->G.Vertex[3].X = x3;
1499 (EleArr + EleCntr - 1)->G.Vertex[3].Y = y3;
1500 (EleArr + EleCntr - 1)->G.Vertex[3].Z = z3;
1501
1502 if (OptElementFiles) {
1503 fprintf(fElem, "##Element Counter: %d\n", EleCntr);
1504 fprintf(fElem, "#DevNb\tCompNb\tPrimNb\tId\n");
1505 fprintf(fElem, "%d\t%d\t%d\t%d\n", (EleArr + EleCntr - 1)->DeviceNb,
1506 (EleArr + EleCntr - 1)->ComponentNb,
1507 (EleArr + EleCntr - 1)->PrimitiveNb,
1508 (EleArr + EleCntr - 1)->Id);
1509 fprintf(fElem, "#GType\tX\tY\tZ\tLX\tLZ\tdA\n");
1510 fprintf(
1511 fElem, "%d\t%.16lg\t%.16lg\t%.16lg\t%.16lg\t%.16lg\t%.16lg\n",
1512 (EleArr + EleCntr - 1)->G.Type, (EleArr + EleCntr - 1)->G.Origin.X,
1513 (EleArr + EleCntr - 1)->G.Origin.Y,
1514 (EleArr + EleCntr - 1)->G.Origin.Z, (EleArr + EleCntr - 1)->G.LX,
1515 (EleArr + EleCntr - 1)->G.LZ, (EleArr + EleCntr - 1)->G.dA);
1516 fprintf(fElem, "#DirnCosn: \n");
1517 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.XUnit.X,
1518 (EleArr + EleCntr - 1)->G.DC.XUnit.Y,
1519 (EleArr + EleCntr - 1)->G.DC.XUnit.Z);
1520 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.YUnit.X,
1521 (EleArr + EleCntr - 1)->G.DC.YUnit.Y,
1522 (EleArr + EleCntr - 1)->G.DC.YUnit.Z);
1523 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.ZUnit.X,
1524 (EleArr + EleCntr - 1)->G.DC.ZUnit.Y,
1525 (EleArr + EleCntr - 1)->G.DC.ZUnit.Z);
1526 fprintf(fElem, "#EType\tLambda\n");
1527 fprintf(fElem, "%d\t%lg\n", (EleArr + EleCntr - 1)->E.Type,
1528 (EleArr + EleCntr - 1)->E.Lambda);
1529 fprintf(fElem, "#NbBCs\tCPX\tCPY\tCPZ\tValue\n");
1530 fprintf(fElem, "%d\t%.16lg\t%.16lg\t%.16lg\t%lg\n",
1531 (EleArr + EleCntr - 1)->BC.NbOfBCs,
1532 (EleArr + EleCntr - 1)->BC.CollPt.X,
1533 (EleArr + EleCntr - 1)->BC.CollPt.Y,
1534 (EleArr + EleCntr - 1)->BC.CollPt.Z,
1535 (EleArr + EleCntr - 1)->BC.Value);
1536 } // if OptElementFiles
1537
1538 // draw centroid and mesh
1539 if (OptGnuplot && OptGnuplotElements) {
1540 fprintf(fgpElem, "%g\t%g\t%g\n", (EleArr + EleCntr - 1)->BC.CollPt.X,
1541 (EleArr + EleCntr - 1)->BC.CollPt.Y,
1542 (EleArr + EleCntr - 1)->BC.CollPt.Z);
1543 } // if OptGnuplot && OptGnuplotElements
1544
1545 if (OptGnuplot && OptGnuplotElements) {
1546 fprintf(fgpMesh, "%g\t%g\t%g\n\n", x0, y0, z0);
1547 fprintf(fgpMesh, "%g\t%g\t%g\n\n", x1, y1, z1);
1548 fprintf(fgpMesh, "%g\t%g\t%g\n\n", x2, y2, z2);
1549 fprintf(fgpMesh, "%g\t%g\t%g\n\n", x3, y3, z3);
1550 fprintf(fgpMesh, "%g\t%g\t%g\n\n", x0, y0, z0);
1551 } // if OptGnuplot && OptGnuplotElements
1552 } // for i
1553 } // for k
1554 ElementEnd[prim] = EleCntr;
1555 NbElmntsOnPrim[prim] = ElementEnd[prim] - ElementBgn[prim] + 1;
1556
1557 if (OptPrimitiveFiles) {
1558 fprintf(fPrim, "Element begin: %d, Element end: %d\n", ElementBgn[prim],
1559 ElementEnd[prim]);
1560 fprintf(fPrim, "Number of elements on primitive: %d\n",
1561 NbElmntsOnPrim[prim]);
1562 fclose(fPrim);
1563 }
1564
1565 if (OptGnuplot && OptGnuplotElements) {
1566 if (prim == 1)
1567 fprintf(fgnuElem, " '%s\' w p", gpElem);
1568 else
1569 fprintf(fgnuElem, ", \\\n \'%s\' w p", gpElem);
1570 if (prim == 1) {
1571 fprintf(fgnuMesh, " '%s\' w l", gpMesh);
1572 fprintf(fgnuMesh, ", \\\n \'%s\' w p ps 1", gpElem);
1573 } else {
1574 fprintf(fgnuMesh, ", \\\n \'%s\' w l", gpMesh);
1575 fprintf(fgnuMesh, ", \\\n \'%s\' w p ps 1", gpElem);
1576 }
1577
1578 fclose(fgpElem);
1579 fclose(fgpMesh);
1580 } // if OptGnuplot && OptGnuplotElements
1581
1582 if (OptElementFiles) fclose(fElem);
1583
1584 return (0);
1585} // end of DiscretizeTriangles
PrintDirnCosn3D
int PrintDirnCosn3D(DirnCosn3D A)
Definition: Vector.c:269
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
DebugLevel
neBEMGLOBAL int DebugLevel
Definition: neBEM.h:235
Epsilon1
neBEMGLOBAL double * Epsilon1
Definition: neBEM.h:75
VolRef1
neBEMGLOBAL int * VolRef1
Definition: neBEM.h:76
VolRef2
neBEMGLOBAL int * VolRef2
Definition: neBEM.h:76
Epsilon2
neBEMGLOBAL double * Epsilon2
Definition: neBEM.h:75

Referenced by SurfaceElements().

◆ DiscretizeWire()

int DiscretizeWire ( 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  NbSegs 
)

Definition at line 409 of file ReTriM.c.

412 {
413 int WireParentObj, WireEType;
414 double WireR, WireL;
415 double WireLambda, WireV;
416 double WireElX, WireElY, WireElZ, WireElL;
417 DirnCosn3D PrimDirnCosn; // direction cosine of the current primitive
418 char primstr[10];
419 char gpElem[256];
420 FILE *fPrim, *fElem, *fgpPrim, *fgpElem;
421
422 // Check inputs
423 if (PrimType[prim] != 2) {
424 neBEMMessage("DiscretizeWire - PrimType in DiscretizeWire");
425 return -1;
426 }
427 if (nvertex != 2) {
428 neBEMMessage("DiscretizeWire - nvertex in DiscretizeWire");
429 return -1;
430 }
431 if (radius < MINDIST) {
432 neBEMMessage("DiscretizeWire - radius in DiscretizeWire");
433 return -1;
434 }
435 if (NbSegs <= 0) {
436 neBEMMessage("DiscretizeWire - NbSegs in DiscretizeWire");
437 return -1;
438 }
439
440 if (OptPrintVertexAndNormal) {
441 printf("nvertex: %d\n", nvertex);
442 for (int vert = 0; vert < nvertex; ++vert) {
443 printf("vert: %d, x: %lg, y: %lg, z: %lg\n", vert, xvert[vert],
444 yvert[vert], zvert[vert]);
445 }
446 printf("radius: %lg\n", radius);
447 } // if OptPrintVertexAndNormal
448
449 // necessary for separating filenames
450 sprintf(primstr, "%d", prim);
451
452 // in order to avoid warning messages
453 fPrim = NULL;
454 fElem = NULL;
455 fgpElem = NULL;
456
457 WireParentObj = 1; // ParentObj not being used now
458
459 WireL = sqrt((xvert[1] - xvert[0]) * (xvert[1] - xvert[0]) // length of wire
460 + (yvert[1] - yvert[0]) * (yvert[1] - yvert[0]) +
461 (zvert[1] - zvert[0]) * (zvert[1] - zvert[0]));
462 WireR = radius;
463 WireElL = WireL / NbSegs; // length of each wire element
464
465 // Direction cosines along the wire - note difference from surface primitives!
466 // The direction along the wire is considered to be the z axis of the LCS
467 // So, let us fix that axial vector first
468 PrimDirnCosn.ZUnit.X = (xvert[1] - xvert[0]) / WireL; // useful
469 PrimDirnCosn.ZUnit.Y = (yvert[1] - yvert[0]) / WireL;
470 PrimDirnCosn.ZUnit.Z = (zvert[1] - zvert[0]) / WireL; // useful
471 // Next, let us find out the coefficients of a plane that passes through the
472 // wire centroid and is normal to the axis of the wire. This is basically the
473 // mid-plane of the cylindrical wire
474 // Any vector OR on the plane normal to the axial direction satisfies
475 // \vec{OR} . \vec{OA} = 0
476 // where O is the wire centroid, A is a point on the axis and R is a point on
477 // the cylindrical surface of the wire. \vec{OA} can be easily replaced by the
478 // axial vector that is equivalent to the vector PrimDirnCosn.ZUnit
479 // The equation of the plane can be shown to be:
480 // XCoef * X + YCoef * Y + ZCoef * Z = Const
481 // double XCoef, YCoef, ZCoef, Const; - not needed any more
482 // double rnorm;
483 // Point3D O, R;
484 // Vector3D OR;
485 // O = CreatePoint3D(WireX, WireY, WireZ);
486 {
487 Vector3D XUnit, YUnit, ZUnit;
488 ZUnit.X = PrimDirnCosn.ZUnit.X;
489 ZUnit.Y = PrimDirnCosn.ZUnit.Y;
490 ZUnit.Z = PrimDirnCosn.ZUnit.Z;
491
492 /* old code - abs instead of fabs??!!
493 XCoef = ZUnit.X;
494 YCoef = ZUnit.Y;
495 ZCoef = ZUnit.Z;
496 double WireX = 0.5 * (xvert[1] + xvert[0]);
497 double WireY = 0.5 * (yvert[1] + yvert[0]);
498 double WireZ = 0.5 * (zvert[1] + zvert[0]);
499 Const = WireX * ZUnit.X + WireY * ZUnit.Y + WireZ * ZUnit.Z;
500 if(abs(XCoef) < 1.0e-12) // X can be anything!
501 {
502 XUnit.X = 1.0;
503 XUnit.Y = 0.0;
504 XUnit.Z = 0.0;
505 YUnit = Vector3DCrossProduct(ZUnit, XUnit);
506 }
507 else
508 {
509 // For a point on the above surface where both Y and Z are zero
510 O = CreatePoint3D(WireX, WireY, WireZ);
511 R = CreatePoint3D(Const, 0, 0);
512 // Create the vector joining O and R; find X and Y unit vectors
513 OR = CreateDistanceVector3D(O,R);
514 XUnit = UnitVector3D(OR);
515 YUnit = Vector3DCrossProduct(ZUnit, XUnit);
516 }
517 old code */
518
519 // replaced following Rob's suggestions (used functions instead of direct
520 // evaluation, although the latter is probably faster)
521 // x-Axis: orthogonal in the 2 largest components.
522 if (fabs(ZUnit.X) >= fabs(ZUnit.Z) && fabs(ZUnit.Y) >= fabs(ZUnit.Z)) {
523 XUnit.X = -ZUnit.Y;
524 XUnit.Y = ZUnit.X;
525 XUnit.Z = 0.0;
526 } else if (fabs(ZUnit.X) >= fabs(ZUnit.Y) &&
527 fabs(ZUnit.Z) >= fabs(ZUnit.Y)) {
528 XUnit.X = -ZUnit.Z;
529 XUnit.Y = 0.0;
530 XUnit.Z = ZUnit.X;
531 } else {
532 XUnit.X = 0.0;
533 XUnit.Y = ZUnit.Z;
534 XUnit.Z = -ZUnit.Y;
535 }
536 XUnit = UnitVector3D(&XUnit);
537
538 // y-Axis: vectorial product of axes 1 and 2.
539 YUnit = Vector3DCrossProduct(&ZUnit, &XUnit);
540 YUnit = UnitVector3D(&YUnit);
541 // end of replacement
542
543 PrimDirnCosn.XUnit.X = XUnit.X;
544 PrimDirnCosn.XUnit.Y = XUnit.Y;
545 PrimDirnCosn.XUnit.Z = XUnit.Z;
546 PrimDirnCosn.YUnit.X = YUnit.X;
547 PrimDirnCosn.YUnit.Y = YUnit.Y;
548 PrimDirnCosn.YUnit.Z = YUnit.Z;
549 } // X and Y direction cosines computed
550
551 // primitive direction cosine assignments
552 PrimDC[prim].XUnit.X = PrimDirnCosn.XUnit.X;
553 PrimDC[prim].XUnit.Y = PrimDirnCosn.XUnit.Y;
554 PrimDC[prim].XUnit.Z = PrimDirnCosn.XUnit.Z;
555 PrimDC[prim].YUnit.X = PrimDirnCosn.YUnit.X;
556 PrimDC[prim].YUnit.Y = PrimDirnCosn.YUnit.Y;
557 PrimDC[prim].YUnit.Z = PrimDirnCosn.YUnit.Z;
558 PrimDC[prim].ZUnit.X = PrimDirnCosn.ZUnit.X;
559 PrimDC[prim].ZUnit.Y = PrimDirnCosn.ZUnit.Y;
560 PrimDC[prim].ZUnit.Z = PrimDirnCosn.ZUnit.Z;
561
562 // primitive origin: also the barycenter for a wire element
563 PrimOriginX[prim] = 0.5 * (xvert[0] + xvert[1]);
564 PrimOriginY[prim] = 0.5 * (yvert[0] + yvert[1]);
565 PrimOriginZ[prim] = 0.5 * (zvert[0] + zvert[1]);
566 PrimLX[prim] = WireR; // radius for wire
567 PrimLZ[prim] = WireL; // length of wire
568
569 WireEType = inttype;
570 WireV = potential;
571 WireLambda = lambda;
572
573 // file output for a primitive
574 if (OptPrimitiveFiles) {
575 char OutPrim[256];
576 strcpy(OutPrim, ModelOutDir);
577 strcat(OutPrim, "/Primitives/Primitive");
578 strcat(OutPrim, primstr);
579 strcat(OutPrim, ".out");
580 fPrim = fopen(OutPrim, "w");
581 if (fPrim == NULL) {
582 neBEMMessage("DiscretizeWire - OutPrim");
583 return -1;
584 }
585 fprintf(fPrim, "#prim: %d, nvertex: %d\n", prim, nvertex);
586 fprintf(fPrim, "Node1: %lg\t%lg\t%lg\n", xvert[0], yvert[0], zvert[0]);
587 fprintf(fPrim, "Node2: %lg\t%lg\t%lg\n", xvert[1], yvert[1], zvert[1]);
588 fprintf(fPrim, "PrimOrigin: %lg\t%lg\t%lg\n", PrimOriginX[prim],
589 PrimOriginY[prim], PrimOriginZ[prim]);
590 fprintf(fPrim, "Primitive lengths: %lg\t%lg\n", PrimLX[prim], PrimLZ[prim]);
591 fprintf(fPrim, "#DirnCosn: \n");
592 fprintf(fPrim, "%lg, %lg, %lg\n", PrimDirnCosn.XUnit.X,
593 PrimDirnCosn.XUnit.Y, PrimDirnCosn.XUnit.Z);
594 fprintf(fPrim, "%lg, %lg, %lg\n", PrimDirnCosn.YUnit.X,
595 PrimDirnCosn.YUnit.Y, PrimDirnCosn.YUnit.Z);
596 fprintf(fPrim, "%lg, %lg, %lg\n", PrimDirnCosn.ZUnit.X,
597 PrimDirnCosn.ZUnit.Y, PrimDirnCosn.ZUnit.Z);
598 fprintf(fPrim, "#volref1: %d, volref2: %d\n", volref1, volref2);
599 fprintf(fPrim, "#NbSegs: %d\n", NbSegs);
600 fprintf(fPrim, "#ParentObj: %d\tEType: %d\n", WireParentObj, WireEType);
601 fprintf(fPrim, "#WireR: %lg\tWireL: %lg\n", WireR, WireL);
602 fprintf(fPrim, "#SurfLambda: %lg\tSurfV: %lg\n", WireLambda, WireV);
603 }
604
605 // necessary for gnuplot
606 if (OptGnuplot && OptGnuplotPrimitives) {
607 char gpPrim[256];
608 strcpy(gpPrim, MeshOutDir);
609 strcat(gpPrim, "/GViewDir/gpPrim");
610 strcat(gpPrim, primstr);
611 strcat(gpPrim, ".out");
612 fgpPrim = fopen(gpPrim, "w");
613 if (fgpPrim == NULL) {
614 neBEMMessage("DiscretizeWire - OutgpPrim");
615 return -1;
616 }
617 fprintf(fgpPrim, "%g\t%g\t%g\n\n", xvert[0], yvert[0], zvert[0]);
618 fprintf(fgpPrim, "%g\t%g\t%g\n\n", xvert[1], yvert[1], zvert[1]);
619 fclose(fgpPrim);
620
621 if (prim == 1)
622 fprintf(fgnuPrim, " '%s\' w l", gpPrim);
623 else
624 fprintf(fgnuPrim, ", \\\n \'%s\' w l", gpPrim);
625 }
626
627 // file outputs for elements on primitive
628 if (OptElementFiles) {
629 char OutElem[256];
630 strcpy(OutElem, MeshOutDir);
631 strcat(OutElem, "/Elements/ElemOnPrim");
632 strcat(OutElem, primstr);
633 strcat(OutElem, ".out");
634 fElem = fopen(OutElem, "w");
635 if (fElem == NULL) {
636 neBEMMessage("DiscretizeWire - OutElem");
637 return -1;
638 }
639 }
640 // gnuplot friendly file outputs for elements on primitive
641 if (OptGnuplot && OptGnuplotElements) {
642 strcpy(gpElem, MeshOutDir);
643 strcat(gpElem, "/GViewDir/gpElemOnPrim");
644 strcat(gpElem, primstr);
645 strcat(gpElem, ".out");
646 fgpElem = fopen(gpElem, "w");
647 if (fgpElem == NULL) {
648 neBEMMessage("DiscretizeWire - OutgpElem");
649 if (fElem) fclose(fElem);
650 return -1;
651 }
652 }
653
654 double xincr = (xvert[1] - xvert[0]) / (double)NbSegs;
655 double yincr = (yvert[1] - yvert[0]) / (double)NbSegs;
656 double zincr = (zvert[1] - zvert[0]) / (double)NbSegs;
657
658 ElementBgn[prim] = EleCntr + 1;
659 double xv0, yv0, zv0, xv1, yv1, zv1;
660 for (int seg = 1; seg <= NbSegs; ++seg) {
661 xv0 = xvert[0] + ((double)seg - 1.0) * xincr;
662 yv0 = yvert[0] + ((double)seg - 1.0) * yincr;
663 zv0 = zvert[0] + ((double)seg - 1.0) * zincr;
664 xv1 = xvert[0] + ((double)seg) * xincr;
665 yv1 = yvert[0] + ((double)seg) * yincr;
666 zv1 = zvert[0] + ((double)seg) * zincr;
667 WireElX = xvert[0] + ((double)seg - 1.0) * xincr + 0.5 * xincr;
668 WireElY = yvert[0] + ((double)seg - 1.0) * yincr + 0.5 * yincr;
669 WireElZ = zvert[0] + ((double)seg - 1.0) * zincr + 0.5 * zincr;
670
671 // Assign element values and write in the file
672 // If element counter exceeds the maximum allowed number of elements, warn!
673 ++EleCntr;
674 if (EleCntr > NbElements) {
675 neBEMMessage("DiscretizeWire - EleCntr more than NbElements!");
676 if (fgpElem) fclose(fgpElem);
677 if (fElem) fclose(fElem);
678 return -1;
679 }
680
681 (EleArr + EleCntr - 1)->DeviceNb =
682 1; // At present, there is only one device
683 (EleArr + EleCntr - 1)->ComponentNb = WireParentObj;
684 (EleArr + EleCntr - 1)->PrimitiveNb = prim;
685 (EleArr + EleCntr - 1)->Id = EleCntr;
686 (EleArr + EleCntr - 1)->G.Type = 2; // linear (wire) here
687 (EleArr + EleCntr - 1)->G.Origin.X = WireElX;
688 (EleArr + EleCntr - 1)->G.Origin.Y = WireElY;
689 (EleArr + EleCntr - 1)->G.Origin.Z = WireElZ;
690 (EleArr + EleCntr - 1)->G.Vertex[0].X = xv0;
691 (EleArr + EleCntr - 1)->G.Vertex[0].Y = yv0;
692 (EleArr + EleCntr - 1)->G.Vertex[0].Z = zv0;
693 (EleArr + EleCntr - 1)->G.Vertex[1].X = xv1;
694 (EleArr + EleCntr - 1)->G.Vertex[1].Y = yv1;
695 (EleArr + EleCntr - 1)->G.Vertex[1].Z = zv1;
696 (EleArr + EleCntr - 1)->G.LX = WireR; // radius of the wire element
697 (EleArr + EleCntr - 1)->G.LZ = WireElL; // wire element length
698 (EleArr + EleCntr - 1)->G.dA = 2.0 * MyPI * (EleArr + EleCntr - 1)->G.LX *
699 (EleArr + EleCntr - 1)->G.LZ;
700 (EleArr + EleCntr - 1)->G.DC.XUnit.X = PrimDirnCosn.XUnit.X;
701 (EleArr + EleCntr - 1)->G.DC.XUnit.Y = PrimDirnCosn.XUnit.Y;
702 (EleArr + EleCntr - 1)->G.DC.XUnit.Z = PrimDirnCosn.XUnit.Z;
703 (EleArr + EleCntr - 1)->G.DC.YUnit.X = PrimDirnCosn.YUnit.X;
704 (EleArr + EleCntr - 1)->G.DC.YUnit.Y = PrimDirnCosn.YUnit.Y;
705 (EleArr + EleCntr - 1)->G.DC.YUnit.Z = PrimDirnCosn.YUnit.Z;
706 (EleArr + EleCntr - 1)->G.DC.ZUnit.X = PrimDirnCosn.ZUnit.X;
707 (EleArr + EleCntr - 1)->G.DC.ZUnit.Y = PrimDirnCosn.ZUnit.Y;
708 (EleArr + EleCntr - 1)->G.DC.ZUnit.Z = PrimDirnCosn.ZUnit.Z;
709 (EleArr + EleCntr - 1)->E.Type = WireEType;
710 (EleArr + EleCntr - 1)->E.Lambda = WireLambda;
711 (EleArr + EleCntr - 1)->Solution = 0.0;
712 (EleArr + EleCntr - 1)->Assigned = charge;
713 (EleArr + EleCntr - 1)->BC.NbOfBCs = 1; // assume one BC per element
714 (EleArr + EleCntr - 1)->BC.CollPt.X =
715 (EleArr + EleCntr - 1)->G.Origin.X; // modify
716 (EleArr + EleCntr - 1)->BC.CollPt.Y =
717 (EleArr + EleCntr - 1)->G.Origin.Y; // to be on
718 (EleArr + EleCntr - 1)->BC.CollPt.Z =
719 (EleArr + EleCntr - 1)->G.Origin.Z; // surface?
720
721 // File operations begin
722 // rfw = fwrite(&Ele, sizeof(Element), 1, fpEle);
723 // printf("Return of fwrite is %d\n", rfw);
724
725 if (OptElementFiles) {
726 fprintf(fElem, "##Element Counter: %d\n", EleCntr);
727 fprintf(fElem, "#DevNb\tCompNb\tPrimNb\tId\n");
728 fprintf(fElem, "%d\t%d\t%d\t%d\n", (EleArr + EleCntr - 1)->DeviceNb,
729 (EleArr + EleCntr - 1)->ComponentNb,
730 (EleArr + EleCntr - 1)->PrimitiveNb, (EleArr + EleCntr - 1)->Id);
731 fprintf(fElem, "#GType\tX\tY\tZ\tLX\tLZ\tdA\n");
732 fprintf(fElem, "%d\t%.16lg\t%.16lg\t%.16lg\t%.16lg\t%.16lg\t%.16lg\n",
733 (EleArr + EleCntr - 1)->G.Type,
734 (EleArr + EleCntr - 1)->G.Origin.X,
735 (EleArr + EleCntr - 1)->G.Origin.Y,
736 (EleArr + EleCntr - 1)->G.Origin.Z, (EleArr + EleCntr - 1)->G.LX,
737 (EleArr + EleCntr - 1)->G.LZ, (EleArr + EleCntr - 1)->G.dA);
738 fprintf(fElem, "#DirnCosn: \n");
739 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.XUnit.X,
740 (EleArr + EleCntr - 1)->G.DC.XUnit.Y,
741 (EleArr + EleCntr - 1)->G.DC.XUnit.Z);
742 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.YUnit.X,
743 (EleArr + EleCntr - 1)->G.DC.YUnit.Y,
744 (EleArr + EleCntr - 1)->G.DC.YUnit.Z);
745 fprintf(fElem, "%lg, %lg, %lg\n", (EleArr + EleCntr - 1)->G.DC.ZUnit.X,
746 (EleArr + EleCntr - 1)->G.DC.ZUnit.Y,
747 (EleArr + EleCntr - 1)->G.DC.ZUnit.Z);
748 fprintf(fElem, "#EType\tLambda\n");
749 fprintf(fElem, "%d\t%lg\n", (EleArr + EleCntr - 1)->E.Type,
750 (EleArr + EleCntr - 1)->E.Lambda);
751 fprintf(fElem, "#NbBCs\tCPX\tCPY\tCPZ\tValue\n");
752 fprintf(fElem, "%d\t%.16lg\t%.16lg\t%.16lg\t%lg\n",
753 (EleArr + EleCntr - 1)->BC.NbOfBCs,
754 (EleArr + EleCntr - 1)->BC.CollPt.X,
755 (EleArr + EleCntr - 1)->BC.CollPt.Y,
756 (EleArr + EleCntr - 1)->BC.CollPt.Z,
757 (EleArr + EleCntr - 1)->BC.Value);
758 } // if OptElementFiles
759 //
760 // mark centroid
761 if (OptGnuplot && OptGnuplotElements) {
762 fprintf(fgpElem, "%g\t%g\t%g\n", (EleArr + EleCntr - 1)->BC.CollPt.X,
763 (EleArr + EleCntr - 1)->BC.CollPt.Y,
764 (EleArr + EleCntr - 1)->BC.CollPt.Z);
765 } // if OptElementFiles
766 // File operations end
767 } // seg loop for wire elements
768 ElementEnd[prim] = EleCntr;
769 NbElmntsOnPrim[prim] = ElementEnd[prim] - ElementBgn[prim] + 1;
770
771 if (OptPrimitiveFiles) {
772 fprintf(fPrim, "Element begin: %d, Element end: %d\n", ElementBgn[prim],
773 ElementEnd[prim]);
774 fprintf(fPrim, "Number of elements on primitive: %d\n",
775 NbElmntsOnPrim[prim]);
776 fclose(fPrim);
777 }
778
779 if (OptElementFiles) fclose(fElem);
780 if (OptGnuplot && OptGnuplotElements) fclose(fgpElem);
781
782 return (0);
783} // end of DiscretizeWire
MyPI
#define MyPI
Definition: ReTriM.c:15
UnitVector3D
Vector3D UnitVector3D(Vector3D *v)
Definition: Vector.c:227

Referenced by WireElements().

◆ SurfaceElements()

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 at line 33 of file ReTriM.c.

36 {
37 // Decide the geometry of this primitive - if it is a rectangle, our job is
38 // greatly simplified. To begin with, we check the number of vertices to
39 // take the decision automatically.
40 // Note that a triangle is the next best bet. All other primitive will have to
41 // be reduced to rectangles (as many as possible) and triangles.
42 // Incidentally, a PrimType (determined in neBEMReadGeom (neBEMInterface.c))
43 // determines whether the primitive is a surface (2D) or a wire (1D),
44 // for a given surface, SurfShape determines whether it is a triangle,
45 // rectangle, or any other polygon besides these two (square is, of course, a
46 // special case of rectangle)
47 int fstatus;
48 switch (nvertex) {
49 case 3: // triangle
50 fstatus = DiscretizeTriangle(prim, nvertex, xvert, yvert, zvert, xnorm,
51 ynorm, znorm, volref1, volref2, inttype,
52 potential, charge, lambda, NbSegX, NbSegZ);
53 // assert(fstatus == 0);
54 if (fstatus != 0) {
55 neBEMMessage("SurfaceElements - DiscretizeTriangle");
56 return -1;
57 }
58 break;
59
60 case 4: // rectangle
61 fstatus = DiscretizeRectangle(prim, nvertex, xvert, yvert, zvert, xnorm,
62 ynorm, znorm, volref1, volref2, inttype,
63 potential, charge, lambda, NbSegX, NbSegZ);
64 // assert(fstatus == 0);
65 if (fstatus != 0) {
66 neBEMMessage("SurfaceElements - DiscretizeRectangle");
67 return -1;
68 }
69 break;
70
71 default:
72 printf("nvertex out of bounds in SurfaceElements ... exiting ...\n");
73 exit(-1);
74 }
75
76 return (0);
77} // end of SurfaceElements
DiscretizeRectangle
int DiscretizeRectangle(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:1590
DiscretizeTriangle
int DiscretizeTriangle(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:788

Referenced by neBEMDiscretize().

◆ WireElements()

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

Definition at line 82 of file ReTriM.c.

85 {
86 int fstatus;
87
88 switch (nvertex) {
89 case 2: // wire
90 fstatus =
91 DiscretizeWire(prim, nvertex, xvert, yvert, zvert, radius, volref1,
92 volref2, inttype, potential, charge, lambda, NbSegs);
93 // assert(fstatus == 0);
94 if (fstatus != 0) {
95 neBEMMessage("WireElements - DiscretizeWire");
96 return -1;
97 }
98 break;
99
100 default:
101 printf("nvertex out of bounds in WireElements ... exiting ...\n");
102 exit(-1);
103 }
104
105 return (0);
106} // end of WireElement
DiscretizeWire
int DiscretizeWire(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 NbSegs)
Definition: ReTriM.c:409

Referenced by neBEMDiscretize().