1 /**********
2 Copyright 1991 Regents of the University of California. All rights reserved.
3 Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
4 Author: 1991 David A. Gates, U. C. Berkeley CAD Group
5 **********/
6
7 #include "ngspice/ngspice.h"
8 #include "ngspice/numglobs.h"
9 #include "ngspice/numconst.h"
10 #include "ngspice/numenum.h"
11 #include "ngspice/twomesh.h"
12 #include "ngspice/twodev.h"
13 #include "ngspice/bool.h"
14 #include "twoddefs.h"
15 #include "twodext.h"
16
17
18
19 static void doMobCoeffs(TWOelem *, int);
20 static void resetEvalFlag(TWOdevice *pDevice);
21
22
23 void
TWObuildMesh(TWOdevice * pDevice,TWOdomain * pDomain,TWOelectrode * pElectrode,TWOmaterial * pMaterial)24 TWObuildMesh(TWOdevice *pDevice, TWOdomain *pDomain,
25 TWOelectrode *pElectrode, TWOmaterial *pMaterial)
26 {
27 int xIndex, yIndex, eIndex, index;
28 int elemType;
29 TWOelem *pElem, *pElem1;
30 TWOnode *pNode, *pNode1, *pNextHNode, *pNextVNode, *pNextDNode;
31 TWOnode ***nodeArray = NULL;
32 TWOedge *pEdge;
33 TWOedge ***edgeArrayH = NULL, ***edgeArrayV = NULL;
34 TWOdomain *pD;
35 TWOelectrode *pE;
36 TWOmaterial *pM;
37 BOOLEAN interiorNode;
38 int poiEqn, numEqn, numElem, numNodes, numEdges;
39 int numXNodes = pDevice->numXNodes;
40 int numYNodes = pDevice->numYNodes;
41 double *xScale = pDevice->xScale;
42 double *yScale = pDevice->yScale;
43 #ifdef NOTDEF
44 FILE *meshFile;
45 #endif
46
47 /* Generate work arrays. */
48 XCALLOC(nodeArray, TWOnode **, 1 + numXNodes);
49 for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
50 XCALLOC(nodeArray[xIndex], TWOnode *, 1 + numYNodes);
51 }
52
53 /* Generate the nodes. */
54 for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
55 for (yIndex = 1; yIndex <= numYNodes; yIndex++) {
56 XCALLOC(pNode, TWOnode, 1);
57 pNode->nodeI = xIndex;
58 pNode->nodeJ = yIndex;
59 pNode->poiEqn = 0;
60 nodeArray[xIndex][yIndex] = pNode;
61 }
62 }
63
64 /* Mark the semiconductor/insulator domains. */
65 if (pDomain == NULL) {
66 fprintf(stderr, "Error: domains not defined for device\n");
67 exit(-1);
68 }
69 for (pD = pDomain; pD != NULL; pD = pD->next) {
70 for (pM = pMaterial; pM != NULL; pM = pM->next) {
71 if (pD->material == pM->id) {
72 break;
73 }
74 }
75 elemType = pM->type;
76 for (xIndex = pD->ixLo; xIndex <= pD->ixHi; xIndex++) {
77 for (yIndex = pD->iyLo; yIndex <= pD->iyHi; yIndex++) {
78 pNode = nodeArray[xIndex][yIndex];
79 pNode->nodeType = elemType;
80 }
81 }
82 }
83 /* Now mark all the metallic domains */
84 for (pE = pElectrode; pE != NULL; pE = pE->next) {
85 for (xIndex = pE->ixLo; xIndex <= pE->ixHi; xIndex++) {
86 for (yIndex = pE->iyLo; yIndex <= pE->iyHi; yIndex++) {
87 pNode = nodeArray[xIndex][yIndex];
88 pNode->nodeType = CONTACT;
89 }
90 }
91 }
92 /*
93 * Now mark as NULL any node in the interior of an electrode, i.e. none of
94 * its neighbors is a different material.
95 */
96 for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
97 for (yIndex = 1; yIndex <= numYNodes; yIndex++) {
98 pNode = nodeArray[xIndex][yIndex];
99 if (pNode->nodeType == CONTACT) {
100 interiorNode = TRUE;
101 if (xIndex > 1) {
102 pNode1 = nodeArray[xIndex - 1][yIndex];
103 if (pNode1->nodeType != 0
104 && pNode1->nodeType != CONTACT) {
105 interiorNode = FALSE;
106 }
107 }
108 if (interiorNode && xIndex < numXNodes) {
109 pNode1 = nodeArray[xIndex + 1][yIndex];
110 if (pNode1->nodeType != 0
111 && pNode1->nodeType != CONTACT) {
112 interiorNode = FALSE;
113 }
114 }
115 if (interiorNode && yIndex > 1) {
116 pNode1 = nodeArray[xIndex][yIndex - 1];
117 if (pNode1->nodeType != 0
118 && pNode1->nodeType != CONTACT) {
119 interiorNode = FALSE;
120 }
121 }
122 if (interiorNode && yIndex < numYNodes) {
123 pNode1 = nodeArray[xIndex][yIndex + 1];
124 if (pNode1->nodeType != 0
125 && pNode1->nodeType != CONTACT) {
126 interiorNode = FALSE;
127 }
128 }
129 if (interiorNode) {
130 pNode->nodeType = 0;
131 }
132 }
133 }
134 }
135
136 /* Delete nodes that do not belong to any domain. */
137 numNodes = 0;
138 for (yIndex = 1; yIndex <= numYNodes; yIndex++) {
139 for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
140 pNode = nodeArray[xIndex][yIndex];
141 if (pNode->nodeType == 0) {
142 /* This node doesn't belong to a domain so delete it. */
143 nodeArray[xIndex][yIndex] = NULL;
144 FREE(pNode);
145 } else {
146 numNodes++;
147 }
148 }
149 }
150 pDevice->numNodes = numNodes;
151
152 /* Now relabel any remaining nodes that are part of electrodes. */
153 setupContacts(pDevice, pElectrode, nodeArray);
154
155 /* Done with the nodes. Now construct the elements and the edges. */
156 numEdges = 0;
157
158 /* Generate the horizontal edges and store in a work array. */
159 XCALLOC(edgeArrayH, TWOedge **, numXNodes);
160 for (xIndex = 1; xIndex < numXNodes; xIndex++) {
161 XCALLOC(edgeArrayH[xIndex], TWOedge *, 1 + numYNodes);
162 }
163 for (yIndex = 1; yIndex <= numYNodes; yIndex++) {
164 for (xIndex = 1; xIndex < numXNodes; xIndex++) {
165 pNode = nodeArray[xIndex][yIndex];
166 pNextHNode = nodeArray[xIndex + 1][yIndex];
167 if (pNode != NULL &&
168 pNextHNode != NULL) {
169 XCALLOC(pEdge, TWOedge, 1);
170 numEdges++;
171 edgeArrayH[xIndex][yIndex] = pEdge;
172 }
173 }
174 }
175
176 /* Generate the vertical edges and store in a work array. */
177 XCALLOC(edgeArrayV, TWOedge **, 1 + numXNodes);
178 for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
179 XCALLOC(edgeArrayV[xIndex], TWOedge *, numYNodes);
180 }
181 for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
182 for (yIndex = 1; yIndex < numYNodes; yIndex++) {
183 pNode = nodeArray[xIndex][yIndex];
184 pNextVNode = nodeArray[xIndex][yIndex + 1];
185 if (pNode != NULL &&
186 pNextVNode != NULL) {
187 XCALLOC(pEdge, TWOedge, 1);
188 numEdges++;
189 edgeArrayV[xIndex][yIndex] = pEdge;
190 }
191 }
192 }
193 pDevice->numEdges = numEdges;
194
195 /* Now construct and count the elements and store the node/edge pointers. */
196 numElem = 1;
197 for (yIndex = 1; yIndex < numYNodes; yIndex++) {
198 for (xIndex = 1; xIndex < numXNodes; xIndex++) {
199 pNode = nodeArray[xIndex][yIndex];
200 pNextHNode = nodeArray[xIndex + 1][yIndex];
201 pNextVNode = nodeArray[xIndex][yIndex + 1];
202 pNextDNode = nodeArray[xIndex + 1][yIndex + 1];
203 if (pNode != NULL &&
204 pNextHNode != NULL &&
205 pNextVNode != NULL &&
206 pNextDNode != NULL) {
207 numElem++;
208 XCALLOC(pElem, TWOelem, 1);
209 pElem->pTLNode = pNode;
210 pElem->pTRNode = pNextHNode;
211 pElem->pBLNode = pNextVNode;
212 pElem->pBRNode = pNextDNode;
213 pElem->pTopEdge = edgeArrayH[xIndex][yIndex];
214 pElem->pBotEdge = edgeArrayH[xIndex][yIndex + 1];
215 pElem->pLeftEdge = edgeArrayV[xIndex][yIndex];
216 pElem->pRightEdge = edgeArrayV[xIndex + 1][yIndex];
217 pDevice->elemArray[xIndex][yIndex] = pElem;
218 } else {
219 pDevice->elemArray[xIndex][yIndex] = NULL;
220 }
221 }
222 }
223
224 /* Create and pack the 1D element array. */
225 pDevice->numElems = numElem - 1;
226 XCALLOC(pDevice->elements, TWOelem *, 1 + numElem);
227 numElem = 1;
228 for (yIndex = 1; yIndex < numYNodes; yIndex++) {
229 for (xIndex = 1; xIndex < numXNodes; xIndex++) {
230 pElem = pDevice->elemArray[xIndex][yIndex];
231 if (pElem != NULL) {
232 pDevice->elements[numElem++] = pElem;
233 }
234 }
235 }
236
237 /* Now create back links from elements to nodes. */
238 for (yIndex = 1; yIndex < numYNodes; yIndex++) {
239 for (xIndex = 1; xIndex < numXNodes; xIndex++) {
240 pElem = pDevice->elemArray[xIndex][yIndex];
241 if (pElem != NULL) {
242 pElem->pTLNode->pBRElem = pElem;
243 pElem->pTRNode->pBLElem = pElem;
244 pElem->pBLNode->pTRElem = pElem;
245 pElem->pBRNode->pTLElem = pElem;
246 if (xIndex > 1) {
247 pElem->pLeftElem = pDevice->elemArray[xIndex-1][yIndex];
248 }
249 if (xIndex < numXNodes - 1) {
250 pElem->pRightElem = pDevice->elemArray[xIndex+1][yIndex];
251 }
252 if (yIndex > 1) {
253 pElem->pTopElem = pDevice->elemArray[xIndex][yIndex-1];
254 }
255 if (yIndex < numYNodes - 1) {
256 pElem->pBotElem = pDevice->elemArray[xIndex][yIndex+1];
257 }
258 }
259 }
260 }
261
262 /* Establish the element types using domain info. */
263 for (pD = pDomain; pD != NULL; pD = pD->next) {
264 for (pM = pMaterial; pM != NULL; pM = pM->next) {
265 if (pD->material == pM->id) {
266 break;
267 }
268 }
269 elemType = pM->type;
270 for (yIndex = pD->iyLo; yIndex < pD->iyHi; yIndex++) {
271 for (xIndex = pD->ixLo; xIndex < pD->ixHi; xIndex++) {
272 pElem = pDevice->elemArray[xIndex][yIndex];
273 if (pElem != NULL) {
274 pElem->domain = pD->id;
275 pElem->elemType = elemType;
276 pElem->matlInfo = pM;
277 }
278 }
279 }
280 }
281
282 /* Establish the edge types. */
283 for (yIndex = 1; yIndex < numYNodes; yIndex++) {
284 for (xIndex = 1; xIndex < numXNodes; xIndex++) {
285 pElem = pDevice->elemArray[xIndex][yIndex];
286 if (pElem != NULL) {
287 for (index = 0; index <= 3; index++) {
288 pEdge = pElem->pEdges[index];
289 pNode = pElem->pNodes[index];
290 pNode1 = pElem->pNodes[(index + 1) % 4];
291 pElem1 = pNode1->pElems[index];
292
293 if (pNode->nodeType == CONTACT &&
294 pNode1->nodeType == CONTACT) {
295 /* Contact edge */
296 pEdge->edgeType = CONTACT;
297 } else if (pNode->nodeType == SCHOTTKY &&
298 pNode1->nodeType == SCHOTTKY) {
299 /* Schottky edge */
300 pEdge->edgeType = SCHOTTKY;
301 } else if (pElem1 == NULL) {
302 /* Neumann edge */
303 pEdge->edgeType = pElem->elemType;
304 } else if (pElem->elemType != pElem1->elemType) {
305 /* Interface edge */
306 pEdge->edgeType = INTERFACE;
307 } else { /* Ignore heterojnxns for now */
308 /* Heterojunction or Homojunction edge */
309 pEdge->edgeType = pElem->elemType;
310 }
311 }
312 }
313 }
314 }
315
316 resetEvalFlag(pDevice);
317 /* Set evaluation flags on nodes and edges. */
318 /* Assign the dx and dy terms in the elements while we're at it. */
319 for (yIndex = 1; yIndex < numYNodes; yIndex++) {
320 for (xIndex = 1; xIndex < numXNodes; xIndex++) {
321 pElem = pDevice->elemArray[xIndex][yIndex];
322 if (pElem != NULL) {
323 pElem->dx = xScale[xIndex + 1] - xScale[xIndex];
324 pElem->dy = yScale[yIndex + 1] - yScale[yIndex];
325 pElem->dxOverDy = pElem->dx / pElem->dy;
326 pElem->dyOverDx = pElem->dy / pElem->dx;
327
328 /*
329 * Semiconductor elements take precedence over Insulator elements, so
330 * set them up first.
331 */
332 for (index = 0; index <= 3; index++) {
333 if (pElem->elemType == SEMICON) {
334 pNode = pElem->pNodes[index];
335 if (!pNode->evaluated) {
336 pNode->evaluated = TRUE;
337 pElem->evalNodes[index] = TRUE;
338 } else {
339 pElem->evalNodes[index] = FALSE;
340 }
341 pEdge = pElem->pEdges[index];
342 if (!pEdge->evaluated) {
343 pEdge->evaluated = TRUE;
344 pElem->evalEdges[index] = TRUE;
345 } else {
346 pElem->evalEdges[index] = FALSE;
347 }
348 }
349 }
350 }
351 }
352 }
353
354 /* Do a second setup pass for Insulator elements */
355 /* Do mobility coefficients now, because we set up dx and dy
356 * in the previous pass
357 */
358 for (yIndex = 1; yIndex < numYNodes; yIndex++) {
359 for (xIndex = 1; xIndex < numXNodes; xIndex++) {
360 pElem = pDevice->elemArray[xIndex][yIndex];
361 if (pElem != NULL) {
362 pElem->direction = 0;
363 pElem->channel = 0;
364 pElem->surface = FALSE;
365 for (index = 0; index <= 3; index++) {
366 if (pElem->elemType == SEMICON) {
367 doMobCoeffs( pElem, index );
368 } else if (pElem->elemType == INSULATOR) {
369 pNode = pElem->pNodes[index];
370 if (!pNode->evaluated) {
371 pNode->evaluated = TRUE;
372 pElem->evalNodes[index] = TRUE;
373 } else {
374 pElem->evalNodes[index] = FALSE;
375 }
376 pEdge = pElem->pEdges[index];
377 if (!pEdge->evaluated) {
378 pEdge->evaluated = TRUE;
379 pElem->evalEdges[index] = TRUE;
380 } else {
381 pElem->evalEdges[index] = FALSE;
382 }
383 }
384 }
385 }
386 }
387 }
388
389 /* Set up the equation numbers for nodes. */
390 poiEqn = numEqn = 1;
391 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
392 pElem = pDevice->elements[eIndex];
393 for (index = 0; index <= 3; index++) {
394 if (pElem->evalNodes[index]) {
395 pNode = pElem->pNodes[index];
396 if (pNode->nodeType != CONTACT) {
397 /* First assign potential equation numbers */
398 if (pNode->nodeType != SCHOTTKY) {
399 pNode->poiEqn = poiEqn++;
400 pNode->psiEqn = numEqn++;
401 }
402 /* Now assign carrier-concentration equation numbers */
403 if (pElem->elemType == INSULATOR) {
404 pNode->nEqn = 0;
405 pNode->pEqn = 0;
406 } else {
407 if (OneCarrier) {
408 /* n and p get same number */
409 pNode->nEqn = numEqn;
410 pNode->pEqn = numEqn++;
411 } else {
412 pNode->nEqn = numEqn++;
413 pNode->pEqn = numEqn++;
414 }
415 }
416 } else { /* This is a contact node. */
417 pNode->poiEqn = 0;
418 pNode->psiEqn = 0;
419 pNode->nEqn = 0;
420 pNode->pEqn = 0;
421 }
422 }
423 }
424 }
425 pDevice->dimEquil = poiEqn;
426 pDevice->dimBias = numEqn;
427
428 /* Open and Print Mesh Output File for Debugging */
429 /* Nuked from release version */
430 #ifdef NOTDEF
431 if (!(meshFile = fopen("mesh.out", "w"))) {
432 perror("mesh.out");
433 exit(-1);
434 }
435 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
436 pElem = pDevice->elements[eIndex];
437 for (index = 0; index <= 3; index++) {
438 if (pElem->evalNodes[index]) {
439 pNode = pElem->pNodes[index];
440 fprintf(meshFile, "node: %5d %5d %5d %5d\n", pNode->nodeI,
441 pNode->nodeJ, pNode->poiEqn, pNode->psiEqn);
442 }
443 }
444 }
445 fflush(meshFile);
446 fclose(meshFile);
447 #endif
448
449 /* Delete work arrays. */
450 for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
451 FREE(nodeArray[xIndex]);
452 FREE(edgeArrayV[xIndex]);
453 }
454 for (xIndex = 1; xIndex < numXNodes; xIndex++) {
455 FREE(edgeArrayH[xIndex]);
456 }
457 FREE(nodeArray);
458 FREE(edgeArrayV);
459 FREE(edgeArrayH);
460
461 /*
462 * TWOprnMesh( pDevice );
463 */
464 }
465
466 void
TWOprnMesh(TWOdevice * pDevice)467 TWOprnMesh(TWOdevice *pDevice)
468 {
469 int eIndex, index;
470 TWOelem *pElem;
471 TWOnode *pNode;
472 TWOedge *pEdge;
473 char *name;
474
475 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
476 pElem = pDevice->elements[eIndex];
477 fprintf(stderr, "elem %5d:\n", eIndex);
478 for (index = 0; index <= 3; index++) {
479 if (pElem->evalNodes[index]) {
480 pNode = pElem->pNodes[index];
481 switch (pNode->nodeType) {
482 case SEMICON:
483 name = "semiconductor";
484 break;
485 case INSULATOR:
486 name = "insulator";
487 break;
488 case CONTACT:
489 name = "contact";
490 break;
491 case SCHOTTKY:
492 name = "schottky";
493 break;
494 case INTERFACE:
495 name = "interface";
496 break;
497 default:
498 name = "unknown";
499 break;
500 }
501 fprintf(stderr, "node %5d: %s %5d %5d\n", index, name,
502 pNode->nodeI, pNode->nodeJ);
503 }
504 if (pElem->evalEdges[index]) {
505 pEdge = pElem->pEdges[index];
506 switch (pEdge->edgeType) {
507 case SEMICON:
508 name = "semiconductor";
509 break;
510 case INSULATOR:
511 name = "insulator";
512 break;
513 case CONTACT:
514 name = "contact";
515 break;
516 case SCHOTTKY:
517 name = "schottky";
518 break;
519 case INTERFACE:
520 name = "interface";
521 break;
522 default:
523 name = "unknown";
524 break;
525 }
526 fprintf(stderr, "edge %5d: %s\n", index, name);
527 }
528 }
529 }
530 }
531
532 /*
533 * We have a separate function for this, so that the setup routines can
534 * reset the state pointers without rebuilding the entire mesh.
535 */
536 void
TWOgetStatePointers(TWOdevice * pDevice,int * numStates)537 TWOgetStatePointers(TWOdevice *pDevice, int *numStates)
538 {
539 int eIndex, index;
540 TWOelem *pElem;
541 TWOnode *pNode;
542 TWOedge *pEdge;
543
544 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
545 pElem = pDevice->elements[eIndex];
546 for (index = 0; index <= 3; index++) {
547 if (pElem->evalNodes[index]) {
548 pNode = pElem->pNodes[index];
549 pNode->nodeState = *numStates;
550 *numStates += TWOnumNodeStates;
551 }
552 if (pElem->evalEdges[index]) {
553 pEdge = pElem->pEdges[index];
554 pEdge->edgeState = *numStates;
555 *numStates += TWOnumEdgeStates;
556 }
557 }
558 }
559 }
560
561 /*
562 * This function computes the percentages of the total semiconductor
563 * width of an edge on the negative and positive sides of the edge
564 */
565 static void
doMobCoeffs(TWOelem * pElem,int index)566 doMobCoeffs(TWOelem *pElem, int index)
567 {
568 TWOelem *pNElem;
569 TWOedge *pEdge;
570 double dl1 = 0.0, dl2 = 0.0;
571
572 pNElem = pElem->pElems[ index ];
573 pEdge = pElem->pEdges[ index ];
574
575 /* If neighbor is not SEMICON, assign and return */
576 if ( (pNElem == NULL) || (pNElem->elemType == INSULATOR) ) {
577 if ( index == 0 || index == 3 ) {
578 pEdge->kNeg = 0.0;
579 pEdge->kPos = 1.0;
580 } else {
581 pEdge->kNeg = 1.0;
582 pEdge->kPos = 0.0;
583 }
584 return;
585 }
586
587 /* Find appropriate dimensions of the elements */
588 switch ( index ) {
589 case 0:
590 dl1 = pNElem->dy;
591 dl2 = pElem->dy;
592 break;
593 case 1:
594 dl1 = pElem->dx;
595 dl2 = pNElem->dx;
596 break;
597 case 2:
598 dl1 = pElem->dy;
599 dl2 = pNElem->dy;
600 break;
601 case 3:
602 dl1 = pNElem->dx;
603 dl2 = pElem->dx;
604 break;
605 }
606
607 /* Assign coefficients */
608 pEdge->kNeg = dl1 / (dl1 + dl2);
609 pEdge->kPos = dl2 / (dl1 + dl2);
610
611 }
612
613 static void
resetEvalFlag(TWOdevice * pDevice)614 resetEvalFlag(TWOdevice *pDevice)
615 {
616 int index, eIndex;
617 TWOelem *pElem;
618
619 for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
620 pElem = pDevice->elements[eIndex];
621 for (index = 0; index <= 3; index++) {
622 pElem->pNodes[index]->evaluated = FALSE;
623 pElem->pEdges[index]->evaluated = FALSE;
624 }
625 }
626 }
627