ngspice/src/ciderlib/twod/twomesh.c

627 lines
16 KiB
C

/**********
Copyright 1991 Regents of the University of California. All rights reserved.
Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
Author: 1991 David A. Gates, U. C. Berkeley CAD Group
**********/
#include "ngspice/ngspice.h"
#include "ngspice/numglobs.h"
#include "ngspice/numconst.h"
#include "ngspice/numenum.h"
#include "ngspice/twomesh.h"
#include "ngspice/twodev.h"
#include "ngspice/bool.h"
#include "twoddefs.h"
#include "twodext.h"
static void doMobCoeffs(TWOelem *, int);
static void resetEvalFlag(TWOdevice *pDevice);
void
TWObuildMesh(TWOdevice *pDevice, TWOdomain *pDomain,
TWOelectrode *pElectrode, TWOmaterial *pMaterial)
{
int xIndex, yIndex, eIndex, index;
int elemType;
TWOelem *pElem, *pElem1;
TWOnode *pNode, *pNode1, *pNextHNode, *pNextVNode, *pNextDNode;
TWOnode ***nodeArray = NULL;
TWOedge *pEdge;
TWOedge ***edgeArrayH = NULL, ***edgeArrayV = NULL;
TWOdomain *pD;
TWOelectrode *pE;
TWOmaterial *pM;
BOOLEAN interiorNode;
int poiEqn, numEqn, numElem, numNodes, numEdges;
int numXNodes = pDevice->numXNodes;
int numYNodes = pDevice->numYNodes;
double *xScale = pDevice->xScale;
double *yScale = pDevice->yScale;
#ifdef NOTDEF
FILE *meshFile;
#endif
/* Generate work arrays. */
XCALLOC(nodeArray, TWOnode **, 1 + numXNodes);
for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
XCALLOC(nodeArray[xIndex], TWOnode *, 1 + numYNodes);
}
/* Generate the nodes. */
for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
for (yIndex = 1; yIndex <= numYNodes; yIndex++) {
XCALLOC(pNode, TWOnode, 1);
pNode->nodeI = xIndex;
pNode->nodeJ = yIndex;
pNode->poiEqn = 0;
nodeArray[xIndex][yIndex] = pNode;
}
}
/* Mark the semiconductor/insulator domains. */
if (pDomain == NULL) {
fprintf(stderr, "Error: domains not defined for device\n");
exit(-1);
}
for (pD = pDomain; pD != NULL; pD = pD->next) {
for (pM = pMaterial; pM != NULL; pM = pM->next) {
if (pD->material == pM->id) {
break;
}
}
elemType = pM->type;
for (xIndex = pD->ixLo; xIndex <= pD->ixHi; xIndex++) {
for (yIndex = pD->iyLo; yIndex <= pD->iyHi; yIndex++) {
pNode = nodeArray[xIndex][yIndex];
pNode->nodeType = elemType;
}
}
}
/* Now mark all the metallic domains */
for (pE = pElectrode; pE != NULL; pE = pE->next) {
for (xIndex = pE->ixLo; xIndex <= pE->ixHi; xIndex++) {
for (yIndex = pE->iyLo; yIndex <= pE->iyHi; yIndex++) {
pNode = nodeArray[xIndex][yIndex];
pNode->nodeType = CONTACT;
}
}
}
/*
* Now mark as NULL any node in the interior of an electrode, i.e. none of
* its neighbors is a different material.
*/
for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
for (yIndex = 1; yIndex <= numYNodes; yIndex++) {
pNode = nodeArray[xIndex][yIndex];
if (pNode->nodeType == CONTACT) {
interiorNode = TRUE;
if (xIndex > 1) {
pNode1 = nodeArray[xIndex - 1][yIndex];
if (pNode1->nodeType != 0
&& pNode1->nodeType != CONTACT) {
interiorNode = FALSE;
}
}
if (interiorNode && xIndex < numXNodes) {
pNode1 = nodeArray[xIndex + 1][yIndex];
if (pNode1->nodeType != 0
&& pNode1->nodeType != CONTACT) {
interiorNode = FALSE;
}
}
if (interiorNode && yIndex > 1) {
pNode1 = nodeArray[xIndex][yIndex - 1];
if (pNode1->nodeType != 0
&& pNode1->nodeType != CONTACT) {
interiorNode = FALSE;
}
}
if (interiorNode && yIndex < numYNodes) {
pNode1 = nodeArray[xIndex][yIndex + 1];
if (pNode1->nodeType != 0
&& pNode1->nodeType != CONTACT) {
interiorNode = FALSE;
}
}
if (interiorNode) {
pNode->nodeType = 0;
}
}
}
}
/* Delete nodes that do not belong to any domain. */
numNodes = 0;
for (yIndex = 1; yIndex <= numYNodes; yIndex++) {
for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
pNode = nodeArray[xIndex][yIndex];
if (pNode->nodeType == 0) {
/* This node doesn't belong to a domain so delete it. */
nodeArray[xIndex][yIndex] = NULL;
FREE(pNode);
} else {
numNodes++;
}
}
}
pDevice->numNodes = numNodes;
/* Now relabel any remaining nodes that are part of electrodes. */
setupContacts(pDevice, pElectrode, nodeArray);
/* Done with the nodes. Now construct the elements and the edges. */
numEdges = 0;
/* Generate the horizontal edges and store in a work array. */
XCALLOC(edgeArrayH, TWOedge **, numXNodes);
for (xIndex = 1; xIndex < numXNodes; xIndex++) {
XCALLOC(edgeArrayH[xIndex], TWOedge *, 1 + numYNodes);
}
for (yIndex = 1; yIndex <= numYNodes; yIndex++) {
for (xIndex = 1; xIndex < numXNodes; xIndex++) {
pNode = nodeArray[xIndex][yIndex];
pNextHNode = nodeArray[xIndex + 1][yIndex];
if (pNode != NULL &&
pNextHNode != NULL) {
XCALLOC(pEdge, TWOedge, 1);
numEdges++;
edgeArrayH[xIndex][yIndex] = pEdge;
}
}
}
/* Generate the vertical edges and store in a work array. */
XCALLOC(edgeArrayV, TWOedge **, 1 + numXNodes);
for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
XCALLOC(edgeArrayV[xIndex], TWOedge *, numYNodes);
}
for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
for (yIndex = 1; yIndex < numYNodes; yIndex++) {
pNode = nodeArray[xIndex][yIndex];
pNextVNode = nodeArray[xIndex][yIndex + 1];
if (pNode != NULL &&
pNextVNode != NULL) {
XCALLOC(pEdge, TWOedge, 1);
numEdges++;
edgeArrayV[xIndex][yIndex] = pEdge;
}
}
}
pDevice->numEdges = numEdges;
/* Now construct and count the elements and store the node/edge pointers. */
numElem = 1;
for (yIndex = 1; yIndex < numYNodes; yIndex++) {
for (xIndex = 1; xIndex < numXNodes; xIndex++) {
pNode = nodeArray[xIndex][yIndex];
pNextHNode = nodeArray[xIndex + 1][yIndex];
pNextVNode = nodeArray[xIndex][yIndex + 1];
pNextDNode = nodeArray[xIndex + 1][yIndex + 1];
if (pNode != NULL &&
pNextHNode != NULL &&
pNextVNode != NULL &&
pNextDNode != NULL) {
numElem++;
XCALLOC(pElem, TWOelem, 1);
pElem->pTLNode = pNode;
pElem->pTRNode = pNextHNode;
pElem->pBLNode = pNextVNode;
pElem->pBRNode = pNextDNode;
pElem->pTopEdge = edgeArrayH[xIndex][yIndex];
pElem->pBotEdge = edgeArrayH[xIndex][yIndex + 1];
pElem->pLeftEdge = edgeArrayV[xIndex][yIndex];
pElem->pRightEdge = edgeArrayV[xIndex + 1][yIndex];
pDevice->elemArray[xIndex][yIndex] = pElem;
} else {
pDevice->elemArray[xIndex][yIndex] = NULL;
}
}
}
/* Create and pack the 1D element array. */
pDevice->numElems = numElem - 1;
XCALLOC(pDevice->elements, TWOelem *, 1 + numElem);
numElem = 1;
for (yIndex = 1; yIndex < numYNodes; yIndex++) {
for (xIndex = 1; xIndex < numXNodes; xIndex++) {
pElem = pDevice->elemArray[xIndex][yIndex];
if (pElem != NULL) {
pDevice->elements[numElem++] = pElem;
}
}
}
/* Now create back links from elements to nodes. */
for (yIndex = 1; yIndex < numYNodes; yIndex++) {
for (xIndex = 1; xIndex < numXNodes; xIndex++) {
pElem = pDevice->elemArray[xIndex][yIndex];
if (pElem != NULL) {
pElem->pTLNode->pBRElem = pElem;
pElem->pTRNode->pBLElem = pElem;
pElem->pBLNode->pTRElem = pElem;
pElem->pBRNode->pTLElem = pElem;
if (xIndex > 1) {
pElem->pLeftElem = pDevice->elemArray[xIndex-1][yIndex];
}
if (xIndex < numXNodes - 1) {
pElem->pRightElem = pDevice->elemArray[xIndex+1][yIndex];
}
if (yIndex > 1) {
pElem->pTopElem = pDevice->elemArray[xIndex][yIndex-1];
}
if (yIndex < numYNodes - 1) {
pElem->pBotElem = pDevice->elemArray[xIndex][yIndex+1];
}
}
}
}
/* Establish the element types using domain info. */
for (pD = pDomain; pD != NULL; pD = pD->next) {
for (pM = pMaterial; pM != NULL; pM = pM->next) {
if (pD->material == pM->id) {
break;
}
}
elemType = pM->type;
for (yIndex = pD->iyLo; yIndex < pD->iyHi; yIndex++) {
for (xIndex = pD->ixLo; xIndex < pD->ixHi; xIndex++) {
pElem = pDevice->elemArray[xIndex][yIndex];
if (pElem != NULL) {
pElem->domain = pD->id;
pElem->elemType = elemType;
pElem->matlInfo = pM;
}
}
}
}
/* Establish the edge types. */
for (yIndex = 1; yIndex < numYNodes; yIndex++) {
for (xIndex = 1; xIndex < numXNodes; xIndex++) {
pElem = pDevice->elemArray[xIndex][yIndex];
if (pElem != NULL) {
for (index = 0; index <= 3; index++) {
pEdge = pElem->pEdges[index];
pNode = pElem->pNodes[index];
pNode1 = pElem->pNodes[(index + 1) % 4];
pElem1 = pNode1->pElems[index];
if (pNode->nodeType == CONTACT &&
pNode1->nodeType == CONTACT) {
/* Contact edge */
pEdge->edgeType = CONTACT;
} else if (pNode->nodeType == SCHOTTKY &&
pNode1->nodeType == SCHOTTKY) {
/* Schottky edge */
pEdge->edgeType = SCHOTTKY;
} else if (pElem1 == NULL) {
/* Neumann edge */
pEdge->edgeType = pElem->elemType;
} else if (pElem->elemType != pElem1->elemType) {
/* Interface edge */
pEdge->edgeType = INTERFACE;
} else { /* Ignore heterojnxns for now */
/* Heterojunction or Homojunction edge */
pEdge->edgeType = pElem->elemType;
}
}
}
}
}
resetEvalFlag(pDevice);
/* Set evaluation flags on nodes and edges. */
/* Assign the dx and dy terms in the elements while we're at it. */
for (yIndex = 1; yIndex < numYNodes; yIndex++) {
for (xIndex = 1; xIndex < numXNodes; xIndex++) {
pElem = pDevice->elemArray[xIndex][yIndex];
if (pElem != NULL) {
pElem->dx = xScale[xIndex + 1] - xScale[xIndex];
pElem->dy = yScale[yIndex + 1] - yScale[yIndex];
pElem->dxOverDy = pElem->dx / pElem->dy;
pElem->dyOverDx = pElem->dy / pElem->dx;
/*
* Semiconductor elements take precedence over Insulator elements, so
* set them up first.
*/
for (index = 0; index <= 3; index++) {
if (pElem->elemType == SEMICON) {
pNode = pElem->pNodes[index];
if (!pNode->evaluated) {
pNode->evaluated = TRUE;
pElem->evalNodes[index] = TRUE;
} else {
pElem->evalNodes[index] = FALSE;
}
pEdge = pElem->pEdges[index];
if (!pEdge->evaluated) {
pEdge->evaluated = TRUE;
pElem->evalEdges[index] = TRUE;
} else {
pElem->evalEdges[index] = FALSE;
}
}
}
}
}
}
/* Do a second setup pass for Insulator elements */
/* Do mobility coefficients now, because we set up dx and dy
* in the previous pass
*/
for (yIndex = 1; yIndex < numYNodes; yIndex++) {
for (xIndex = 1; xIndex < numXNodes; xIndex++) {
pElem = pDevice->elemArray[xIndex][yIndex];
if (pElem != NULL) {
pElem->direction = 0;
pElem->channel = 0;
pElem->surface = FALSE;
for (index = 0; index <= 3; index++) {
if (pElem->elemType == SEMICON) {
doMobCoeffs( pElem, index );
} else if (pElem->elemType == INSULATOR) {
pNode = pElem->pNodes[index];
if (!pNode->evaluated) {
pNode->evaluated = TRUE;
pElem->evalNodes[index] = TRUE;
} else {
pElem->evalNodes[index] = FALSE;
}
pEdge = pElem->pEdges[index];
if (!pEdge->evaluated) {
pEdge->evaluated = TRUE;
pElem->evalEdges[index] = TRUE;
} else {
pElem->evalEdges[index] = FALSE;
}
}
}
}
}
}
/* Set up the equation numbers for nodes. */
poiEqn = numEqn = 1;
for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
pElem = pDevice->elements[eIndex];
for (index = 0; index <= 3; index++) {
if (pElem->evalNodes[index]) {
pNode = pElem->pNodes[index];
if (pNode->nodeType != CONTACT) {
/* First assign potential equation numbers */
if (pNode->nodeType != SCHOTTKY) {
pNode->poiEqn = poiEqn++;
pNode->psiEqn = numEqn++;
}
/* Now assign carrier-concentration equation numbers */
if (pElem->elemType == INSULATOR) {
pNode->nEqn = 0;
pNode->pEqn = 0;
} else {
if (OneCarrier) {
/* n and p get same number */
pNode->nEqn = numEqn;
pNode->pEqn = numEqn++;
} else {
pNode->nEqn = numEqn++;
pNode->pEqn = numEqn++;
}
}
} else { /* This is a contact node. */
pNode->poiEqn = 0;
pNode->psiEqn = 0;
pNode->nEqn = 0;
pNode->pEqn = 0;
}
}
}
}
pDevice->dimEquil = poiEqn;
pDevice->dimBias = numEqn;
/* Open and Print Mesh Output File for Debugging */
/* Nuked from release version */
#ifdef NOTDEF
if (!(meshFile = fopen("mesh.out", "w"))) {
perror("mesh.out");
exit(-1);
}
for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
pElem = pDevice->elements[eIndex];
for (index = 0; index <= 3; index++) {
if (pElem->evalNodes[index]) {
pNode = pElem->pNodes[index];
fprintf(meshFile, "node: %5d %5d %5d %5d\n", pNode->nodeI,
pNode->nodeJ, pNode->poiEqn, pNode->psiEqn);
}
}
}
fflush(meshFile);
fclose(meshFile);
#endif
/* Delete work arrays. */
for (xIndex = 1; xIndex <= numXNodes; xIndex++) {
FREE(nodeArray[xIndex]);
FREE(edgeArrayV[xIndex]);
}
for (xIndex = 1; xIndex < numXNodes; xIndex++) {
FREE(edgeArrayH[xIndex]);
}
FREE(nodeArray);
FREE(edgeArrayV);
FREE(edgeArrayH);
/*
* TWOprnMesh( pDevice );
*/
}
void
TWOprnMesh(TWOdevice *pDevice)
{
int eIndex, index;
TWOelem *pElem;
TWOnode *pNode;
TWOedge *pEdge;
char *name;
for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
pElem = pDevice->elements[eIndex];
fprintf(stderr, "elem %5d:\n", eIndex);
for (index = 0; index <= 3; index++) {
if (pElem->evalNodes[index]) {
pNode = pElem->pNodes[index];
switch (pNode->nodeType) {
case SEMICON:
name = "semiconductor";
break;
case INSULATOR:
name = "insulator";
break;
case CONTACT:
name = "contact";
break;
case SCHOTTKY:
name = "schottky";
break;
case INTERFACE:
name = "interface";
break;
default:
name = "unknown";
break;
}
fprintf(stderr, "node %5d: %s %5d %5d\n", index, name,
pNode->nodeI, pNode->nodeJ);
}
if (pElem->evalEdges[index]) {
pEdge = pElem->pEdges[index];
switch (pEdge->edgeType) {
case SEMICON:
name = "semiconductor";
break;
case INSULATOR:
name = "insulator";
break;
case CONTACT:
name = "contact";
break;
case SCHOTTKY:
name = "schottky";
break;
case INTERFACE:
name = "interface";
break;
default:
name = "unknown";
break;
}
fprintf(stderr, "edge %5d: %s\n", index, name);
}
}
}
}
/*
* We have a separate function for this, so that the setup routines can
* reset the state pointers without rebuilding the entire mesh.
*/
void
TWOgetStatePointers(TWOdevice *pDevice, int *numStates)
{
int eIndex, index;
TWOelem *pElem;
TWOnode *pNode;
TWOedge *pEdge;
for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
pElem = pDevice->elements[eIndex];
for (index = 0; index <= 3; index++) {
if (pElem->evalNodes[index]) {
pNode = pElem->pNodes[index];
pNode->nodeState = *numStates;
*numStates += TWOnumNodeStates;
}
if (pElem->evalEdges[index]) {
pEdge = pElem->pEdges[index];
pEdge->edgeState = *numStates;
*numStates += TWOnumEdgeStates;
}
}
}
}
/*
* This function computes the percentages of the total semiconductor
* width of an edge on the negative and positive sides of the edge
*/
static void
doMobCoeffs(TWOelem *pElem, int index)
{
TWOelem *pNElem;
TWOedge *pEdge;
double dl1 = 0.0, dl2 = 0.0;
pNElem = pElem->pElems[ index ];
pEdge = pElem->pEdges[ index ];
/* If neighbor is not SEMICON, assign and return */
if ( (pNElem == NULL) || (pNElem->elemType == INSULATOR) ) {
if ( index == 0 || index == 3 ) {
pEdge->kNeg = 0.0;
pEdge->kPos = 1.0;
} else {
pEdge->kNeg = 1.0;
pEdge->kPos = 0.0;
}
return;
}
/* Find appropriate dimensions of the elements */
switch ( index ) {
case 0:
dl1 = pNElem->dy;
dl2 = pElem->dy;
break;
case 1:
dl1 = pElem->dx;
dl2 = pNElem->dx;
break;
case 2:
dl1 = pElem->dy;
dl2 = pNElem->dy;
break;
case 3:
dl1 = pNElem->dx;
dl2 = pElem->dx;
break;
}
/* Assign coefficients */
pEdge->kNeg = dl1 / (dl1 + dl2);
pEdge->kPos = dl2 / (dl1 + dl2);
}
static void
resetEvalFlag(TWOdevice *pDevice)
{
int index, eIndex;
TWOelem *pElem;
for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
pElem = pDevice->elements[eIndex];
for (index = 0; index <= 3; index++) {
pElem->pNodes[index]->evaluated = FALSE;
pElem->pEdges[index]->evaluated = FALSE;
}
}
}