You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
885 lines
25 KiB
885 lines
25 KiB
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
#include "portmesh.h"
|
|
#include "material.h"
|
|
#include "bc.h"
|
|
#include "tensor.h"
|
|
#include "cvtr.h"
|
|
#include "Constants.h"
|
|
#include "tri.elem"
|
|
#include "vtkwriter.h"
|
|
|
|
portMesh::portMesh(){
|
|
faceCNT = 0;
|
|
edgeCNT = 0;
|
|
nodeCNT = 0;
|
|
objCNT = 0;
|
|
objMAP = 0;
|
|
ndArray = 0;
|
|
fcArray = 0;
|
|
modeIMP = 0.;
|
|
magE = 0.;
|
|
impZ = 0.;
|
|
globToLocMap_ = 0;
|
|
}
|
|
|
|
portMesh::~portMesh(){
|
|
delete [] fcArray;
|
|
delete [] ndArray;
|
|
delete [] objMAP;
|
|
delete globToLocMap_;
|
|
}
|
|
|
|
void portMesh::writeMesh(Material *propMAT, fp_t freq)
|
|
{
|
|
FILE *foo;
|
|
char meshName[180];
|
|
int i, j, k, objNum;
|
|
tensor epsr, mur;
|
|
node *nd;
|
|
vtr vv, vv2D;
|
|
|
|
sprintf(meshName, "%s.mesh", portName);
|
|
foo = fopen(meshName, "w");
|
|
fprintf(foo, "%.20e %.20e %.20e\n", (coord.O.getx()), (coord.O.gety()), (coord.O.getz()));
|
|
fprintf(foo, "%.20e %.20e %.20e\n", (coord.x_axis.getx()), (coord.x_axis.gety()), (coord.x_axis.getz()));
|
|
fprintf(foo, "%.20e %.20e %.20e\n", (coord.y_axis.getx()), (coord.y_axis.gety()), (coord.y_axis.getz()));
|
|
fprintf(foo, "1.0 \n");
|
|
fprintf(foo, "%d %d %d \n", objCNT, nodeCNT, faceCNT);
|
|
|
|
for(i = 0; i < nodeCNT; i++)
|
|
{
|
|
nd = ndArray[i];
|
|
vv = nd->getCoord();
|
|
vv2D = coord.Transform(vv);
|
|
|
|
fprintf(foo, "0 %.20e %.20e \n", (vv2D.getx()), (vv2D.gety()));
|
|
}
|
|
|
|
for(i = 0; i < faceCNT; i++){
|
|
objNum = fcArray[i]->gettetraPtr(0)->getobjNum();
|
|
fprintf(foo, "%d %d %d %d \n", getobjMAP(objNum) + 1, fcArray[i]->getNodePtr(0)->num2D, fcArray[i]->getNodePtr(1)->num2D, fcArray[i]->getNodePtr(2)->num2D);
|
|
for(j = 0; j < 3; j++){
|
|
if(fcArray[i]->getEdge(j) == 0)
|
|
fprintf(foo, " 1 ");
|
|
else
|
|
fprintf(foo, " 0 ");
|
|
}
|
|
|
|
fprintf(foo, "\n");
|
|
}
|
|
|
|
for(i = 0; i < objCNT; i++){
|
|
objNum = objMAP[i];
|
|
epsr = propMAT[objNum].epsr;
|
|
mur = propMAT[objNum].mur;
|
|
epsr = coord.Transform(epsr);
|
|
|
|
fprintf(foo, "%d \n", i + 1);
|
|
for(j = 0; j < 3; j++){
|
|
for(k = 0; k < 3; k++){
|
|
fprintf(foo, "%.20e ", (epsr.getEntry(j, k)));
|
|
}
|
|
fprintf(foo, "\n");
|
|
}
|
|
|
|
for(j = 0; j < 3; j++){
|
|
for(k = 0; k < 3; k++){
|
|
fprintf(foo, "%.20e ", (mur.getEntry(j, k)));
|
|
}
|
|
fprintf(foo, "\n");
|
|
}
|
|
|
|
for(j = 0; j < 3; j++){
|
|
for(k = 0; k < 3; k++){
|
|
fprintf(foo, "0.0 ");
|
|
}
|
|
fprintf(foo, "\n");
|
|
}
|
|
|
|
for(j = 0; j < 3; j++){
|
|
for(k = 0; k < 3; k++){
|
|
fprintf(foo, "0.0 ");
|
|
}
|
|
fprintf(foo, "\n");
|
|
}
|
|
}
|
|
|
|
fprintf(foo, "%d \n", vline.np);
|
|
for(i = 0; i < vline.np; i++){
|
|
vv = vline.coord[i];
|
|
cout << vv.getx() << " " << vv.gety() << " " << vv.getz() << endl;
|
|
vv2D = coord.Transform(vv);
|
|
fprintf(foo, "%.20e %.20e \n", vv2D.getx(), vv2D.gety());
|
|
}
|
|
fclose(foo);
|
|
}
|
|
|
|
|
|
void portMesh::writeMesh(Material* propMat){
|
|
int i, j, k;
|
|
vtr vv, vv2D;
|
|
|
|
char meshName[STRLEN];
|
|
memset(meshName, 0, STRLEN*sizeof(char));
|
|
sprintf(meshName, "%s.mesh", portName);
|
|
FILE* foo = fopen(meshName, "w");
|
|
vtr& O = coord.getO();
|
|
vtr& xAxis = coord.getXAxis();
|
|
vtr& yAxis = coord.getYAxis();
|
|
|
|
fprintf(foo, "%.20e %.20e %.20e\n", (O.getx()), (O.gety()), (O.getz()));
|
|
fprintf(foo, "%.20e %.20e %.20e\n", (xAxis.getx()), (xAxis.gety()), (xAxis.getz()));
|
|
fprintf(foo, "%.20e %.20e %.20e\n", (yAxis.getx()), (yAxis.gety()), (yAxis.getz()));
|
|
fprintf(foo, "1.0\n");
|
|
fprintf(foo, "%d %d %d \n", objCNT, nodeCNT, faceCNT);
|
|
|
|
for(i = 0; i < nodeCNT; i++){
|
|
node* nd = ndArray[i];
|
|
vv = nd->getCoord();
|
|
vv2D = coord.Transform(vv);
|
|
fprintf(foo, "0 %.20e %.20e\n", (vv2D.getx()), (vv2D.gety()));
|
|
}
|
|
|
|
int objNum = 0;
|
|
for(i = 0; i < faceCNT; i++){
|
|
objNum = fcArray[i]->getTetra(0).getobjNum();
|
|
fprintf(foo, "%d %d %d %d\n", getobjMAP(objNum) + 1,
|
|
globToLocMap_->find(fcArray[i]->getNode(0)->getid())->second,
|
|
globToLocMap_->find(fcArray[i]->getNode(1)->getid())->second,
|
|
globToLocMap_->find(fcArray[i]->getNode(2)->getid())->second);
|
|
for(j = 0; j < 3; j++){
|
|
if(fcArray[i]->getEdge(j)->getbType() == pecType)
|
|
fprintf(foo, " 1 ");
|
|
else
|
|
fprintf(foo, " 0 ");
|
|
}
|
|
fprintf(foo, "\n");
|
|
}
|
|
|
|
tensor epsr;
|
|
tensor mur;
|
|
for(i = 0; i < objCNT; i++){
|
|
objNum = objMAP[i];
|
|
epsr = propMat[objNum].epsr;
|
|
mur = propMat[objNum].mur;
|
|
epsr = coord.Transform(epsr);
|
|
|
|
fprintf(foo, "%d\n", i + 1);
|
|
for(j = 0; j < 3; j++){
|
|
for(k = 0; k < 3; k++)
|
|
fprintf(foo, "%.20e ", (epsr.getEntry(j, k)));
|
|
fprintf(foo, "\n");
|
|
}
|
|
for(j = 0; j < 3; j++){
|
|
for(k = 0; k < 3; k++)
|
|
fprintf(foo, "%.20e ", (mur.getEntry(j, k)));
|
|
fprintf(foo, "\n");
|
|
}
|
|
for(j = 0; j < 3; j++){
|
|
for(k = 0; k < 3; k++)
|
|
fprintf(foo, "%.20e ", 0.0);
|
|
fprintf(foo, "\n");
|
|
}
|
|
for(j = 0; j < 3; j++){
|
|
for(k = 0; k < 3; k++) fprintf(foo, "0.0 ");
|
|
fprintf(foo, "\n");
|
|
}
|
|
}
|
|
|
|
fprintf(foo, "%d\n", vline.np);
|
|
|
|
for(i = 0; i < vline.np; i++){
|
|
cout << vv.getx() << " " << vv.gety() << " " << vv.getz() << endl;
|
|
vv = vline.coord[i];
|
|
vv2D = coord.Transform(vv);
|
|
fprintf(foo, "%.20e %.20e\n", (vv2D.getx()), (vv2D.gety()));
|
|
}
|
|
|
|
fclose(foo);
|
|
}
|
|
|
|
void portMesh::makeObjMap(){
|
|
int i, tmpCNT, objNum;
|
|
|
|
tmpCNT = 0;
|
|
for(i = 0; i < faceCNT; i++){
|
|
objNum = fcArray[i]->gettetraPtr(0)->getobjNum();
|
|
tmpCNT = (tmpCNT > objNum) ? tmpCNT : objNum;
|
|
}
|
|
tmpCNT++;
|
|
|
|
objMAP = new int[tmpCNT];
|
|
for(i = 0; i < faceCNT; i++){
|
|
objNum = fcArray[i]->gettetraPtr(0)->getobjNum();
|
|
insertOBJ(objNum);
|
|
}
|
|
}
|
|
|
|
int portMesh::getobjMAP(int obj){
|
|
int i;
|
|
for(i = 0; i < objCNT; i++){
|
|
if(objMAP[i] == obj) return i;
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
void portMesh::insertOBJ(int obj){
|
|
int i;
|
|
|
|
for(i = 0; i < objCNT; i++){
|
|
if(objMAP[i] == obj)
|
|
return;
|
|
}
|
|
|
|
objMAP[objCNT] = obj;
|
|
objCNT++;
|
|
}
|
|
|
|
void portMesh::setFaceCnt(int faceCnt){
|
|
faceCNT = faceCnt;
|
|
fcArray = new face*[faceCNT];
|
|
memset(fcArray, 0, faceCNT * sizeof(face*));
|
|
}
|
|
|
|
void portMesh::setNodeCnt(int nodeCnt){
|
|
nodeCNT = nodeCnt;
|
|
ndArray = new node*[nodeCNT];
|
|
memset(ndArray, 0, nodeCNT * sizeof(node*));
|
|
}
|
|
|
|
// Read from portnems.vpath the definition of the path which defines
|
|
// the voltage integration forthe mode. This is necessary to do
|
|
// de-embedding to 50 ohm transmission lines. ifthe file does not exist
|
|
// then the S parameters will be expressed in terms of the intrinsic
|
|
// impedances.
|
|
//
|
|
|
|
void portMesh::readVline(fp_t unit_conv){
|
|
FILE *foo;
|
|
char vlineName[180];
|
|
int i, np;
|
|
fp_t x, y, z;
|
|
double x_,y_,z_;
|
|
|
|
sprintf(vlineName, "%s.vpath", portName);
|
|
if((foo = fopen(vlineName, "r")) == NULL){
|
|
vline.np = 0;
|
|
vline.coord = 0;
|
|
|
|
return;
|
|
}
|
|
|
|
fscanf(foo, "%d", &np);
|
|
vline.np = np;
|
|
vline.coord = new vtr[np];
|
|
for(i = 0; i < np; i++){
|
|
fscanf(foo, "%le", &x_);
|
|
fscanf(foo, "%le", &y_);
|
|
fscanf(foo, "%le", &z_);
|
|
(vline.coord[i]).setvtr((fp_t)x_ * unit_conv, (fp_t)y_ * unit_conv, (fp_t)z_ * unit_conv);
|
|
}
|
|
fclose(foo);
|
|
}
|
|
|
|
void portMesh::readVline(fp_t unit_conv, fp_t& xoff, fp_t& yoff, fp_t& zoff){
|
|
FILE *foo;
|
|
char vlineName[180];
|
|
int i, np;
|
|
fp_t x, y, z;
|
|
double x_,y_,z_;
|
|
|
|
sprintf(vlineName, "%s.vpath", portName);
|
|
if((foo = fopen(vlineName, "r")) == NULL){
|
|
vline.np = 0;
|
|
vline.coord = 0;
|
|
|
|
return;
|
|
}
|
|
|
|
fscanf(foo, "%d", &np);
|
|
vline.np = np;
|
|
vline.coord = new vtr[np];
|
|
for(i = 0; i < np; i++){
|
|
fscanf(foo, "%le", &x_);
|
|
fscanf(foo, "%le", &y_);
|
|
fscanf(foo, "%le", &z_);
|
|
(vline.coord[i]).setvtr((fp_t)x_ * unit_conv + xoff, (fp_t)y_ * unit_conv + yoff, (fp_t)z_ * unit_conv + zoff);
|
|
}
|
|
fclose(foo);
|
|
}
|
|
|
|
// portMesh::makeRHS
|
|
// foreach port, we assemble a right-hand-side to be used as the
|
|
// excitation forthe computations of the entire scattering matrix
|
|
// of the multiport microwave components.
|
|
// void portMesh::makeRHS(fp_t freq, int dim, int *Part, portCounter pC){
|
|
// FILE *foo;
|
|
// char solname[180], cc;
|
|
// int i, j, k, tmpCNT, triMAP[8];
|
|
// vtr tmpVTR, nvtr[2];
|
|
// vtr hvtr[6], jvtr[6];
|
|
// cVtr tvtr;
|
|
// fp_t jval[8], rhs, cval;
|
|
// face *fc;
|
|
// denseMat elemM(8, 8);
|
|
// Material prop;
|
|
// fp_t scaleF;
|
|
|
|
// cval = -Uo;
|
|
// portY.N = dim;
|
|
// portY.entry = new fp_t[dim];
|
|
|
|
// sprintf(solname, "%s.1H1", portName);
|
|
// if((foo = fopen(solname, "rb")) == NULL){
|
|
// printf("Could not open File %s \n", solname);
|
|
// exit(1);
|
|
// }
|
|
// fread(&tmpCNT, sizeof(int ), 1, foo);
|
|
// fread(&cc, sizeof(char ), 1, foo);
|
|
// fread(&cc, sizeof(char ), 1, foo);
|
|
// fread(&tmpCNT, sizeof(int ), 1, foo);
|
|
// fread(&gamma, sizeof(fp_t ), 1, foo);
|
|
// fread(&gamma, sizeof(fp_t ), 1, foo);
|
|
// for(i = 0; i < 3; i++){
|
|
// fread(&tmpVTR, sizeof(vtr ), 1, foo);
|
|
// }
|
|
// fread(&modeIMP, sizeof(fp_t ), 1, foo);
|
|
// // scaleF = modeIMP;
|
|
// // scaleF = scaleF / impZ;
|
|
// // scaleF = sqrt(scaleF);
|
|
// scaleF = 1.0;
|
|
// for(i = 0; i < 8; i++)
|
|
// jval[i] = 0.0;
|
|
// for(i = 0; i < faceCNT; i++){
|
|
// fc = fcArray[i];
|
|
// maketriMAP(fc, Part, triMAP, pC);
|
|
|
|
// for(j = 0; j < 6; j++){
|
|
// fread(&tvtr, sizeof(cVtr ), 1, foo);
|
|
// hvtr[j].setvtr(tvtr.x.real(), tvtr.y.real(), tvtr.z.real());
|
|
// hvtr[j] = hvtr[j] * scaleF;
|
|
|
|
// jvtr[j].setvtr(hvtr[j].gety() * (-1.0), hvtr[j].getx(), 0.0);
|
|
// }
|
|
// fc->P2toH1(jvtr, jval, coord);
|
|
// fc->makeElemMatrix(prop, coord, 8, &elemM);
|
|
// for(j = 0; j < 8; j++){
|
|
// if(triMAP[j] < 0)
|
|
// continue;
|
|
// rhs = 0.0;
|
|
// for(k = 0; k < 8; k++){
|
|
// rhs = rhs + jval[k] * elemM.getEntry(j, k);
|
|
// }
|
|
// portY.addentry(triMAP[j], rhs * cval);
|
|
// }
|
|
// }
|
|
// fclose(foo);
|
|
// }
|
|
|
|
struct dbl_vtr{
|
|
double x,y,z;
|
|
};
|
|
|
|
struct dbl_cvtr{
|
|
complex<double> x,y,z;
|
|
};
|
|
|
|
|
|
|
|
// -----------------------------------------------------------------------------
|
|
// portMesh::makeRHS_E
|
|
// Reads electric-field port samples from "<portName>.1E1", one block per face,
|
|
// 6 complex vector samples per face (second-order triangle).
|
|
//
|
|
// Each sample dbl_cvtr has .x/.y/.z as std::complex<double>.
|
|
// Samples are rotated into global coordinates via coord.x_axis/y_axis/z_axis
|
|
// and scaled by scaleF (currently 1).
|
|
// -----------------------------------------------------------------------------
|
|
void portMesh::makeRHS_E()
|
|
{
|
|
FILE *foo;
|
|
char solname[StrOutput], cc;
|
|
int i, j, tmpCNT, triMAP[8];
|
|
vtr nvtr[2];
|
|
dbl_vtr tmpVTR;
|
|
vtr evtr[NumOfEdges]; // holds 6 E samples for the current face
|
|
dbl_cvtr tvtr;
|
|
fp_t cval;
|
|
face *fc;
|
|
fp_t scaleF;
|
|
|
|
double gamma_,modeIMP_;
|
|
|
|
sprintf(solname, "%s.1E1", portName);
|
|
|
|
cout << "Reading port mode solution from file: " << solname << endl;
|
|
|
|
|
|
if((foo = fopen(solname, "rb")) == NULL)
|
|
{
|
|
printf("Could not open File %s \n", solname);
|
|
exit(1);
|
|
}
|
|
|
|
|
|
fread(&tmpCNT, sizeof(int), 1, foo);
|
|
fread(&cc, sizeof(char), 1, foo);
|
|
fread(&cc, sizeof(char), 1, foo);
|
|
fread(&tmpCNT, sizeof(int), 1, foo);
|
|
fread(&gamma_, sizeof(double), 1, foo);
|
|
fread(&gamma_, sizeof(double), 1, foo);
|
|
for(i = 0; i < NumOfUnitaryVectors; i++)
|
|
{
|
|
fread(&tmpVTR, sizeof(dbl_vtr), 1, foo);
|
|
}
|
|
fread(&modeIMP_, sizeof(double), 1, foo);
|
|
|
|
gamma = (fp_t)gamma_;
|
|
modeIMP = (fp_t)modeIMP_;
|
|
|
|
// --- per-face blocks ---
|
|
scaleF = 1.0;
|
|
for(j = 0; j < NumOfEdges; j++)
|
|
{
|
|
evtr[j].setvtr(0.0, 0.0, 0.0);
|
|
}
|
|
|
|
// read 6 complex-vector samples (x,y,z) for this face
|
|
for(i = 0; i < faceCNT; i++)
|
|
{
|
|
fc = fcArray[i];
|
|
for(j = 0; j < 6; j++)
|
|
{
|
|
fread(&tvtr, sizeof(dbl_cvtr), 1, foo);
|
|
|
|
// rotate from file-local axes to global axes
|
|
evtr[j] = coord.x_axis * (fp_t)(tvtr.x.real()) + coord.y_axis * (fp_t)(tvtr.y.real()) + coord.z_axis * (fp_t)(tvtr.z.real());
|
|
evtr[j] = evtr[j] * scaleF;
|
|
}
|
|
|
|
// convert the 6 P2 samples to H1 coefficients and store on the face
|
|
fc->Set_WgSolverExcitation_E(evtr);
|
|
}
|
|
fclose(foo);
|
|
}
|
|
|
|
|
|
|
|
void portMesh::makeRHS_H()
|
|
{
|
|
FILE *foo;
|
|
char solname[180], cc;
|
|
int i, j, tmpCNT, triMAP[8];
|
|
vtr nvtr[2];
|
|
dbl_vtr tmpVTR;
|
|
vtr hvtr[6];
|
|
dbl_cvtr tvtr;
|
|
fp_t cval;
|
|
face *fc;
|
|
fp_t scaleF;
|
|
double gamma_,modeIMP_;
|
|
|
|
sprintf(solname, "%s.1H1", portName);
|
|
if((foo = fopen(solname, "rb")) == NULL){
|
|
printf("Could not open File %s \n", solname);
|
|
exit(1);
|
|
}
|
|
fread(&tmpCNT, sizeof(int), 1, foo);
|
|
fread(&cc, sizeof(char), 1, foo);
|
|
fread(&cc, sizeof(char), 1, foo);
|
|
fread(&tmpCNT, sizeof(int), 1, foo);
|
|
fread(&gamma_, sizeof(double), 1, foo);
|
|
fread(&gamma_, sizeof(double), 1, foo);
|
|
|
|
for(i = 0; i < 3; i++)
|
|
{
|
|
fread(&tmpVTR, sizeof(dbl_vtr), 1, foo);
|
|
}
|
|
|
|
fread(&modeIMP_, sizeof(double ), 1, foo);
|
|
|
|
gamma = (fp_t)gamma_;
|
|
modeIMP = (fp_t)modeIMP_;
|
|
|
|
scaleF = 1.0;
|
|
for(j = 0; j < 6; j++){
|
|
hvtr[j].setvtr(0.0, 0.0, 0.0);
|
|
}
|
|
for(i = 0; i < faceCNT; i++){
|
|
fc = fcArray[i];
|
|
for(j = 0; j < 6; j++){
|
|
fread(&tvtr, sizeof(dbl_cvtr), 1, foo);
|
|
hvtr[j] = coord.x_axis * (fp_t)(tvtr.x.real()) + coord.y_axis * (fp_t)(tvtr.y.real()) + coord.z_axis * (fp_t)(tvtr.z.real());
|
|
hvtr[j] = hvtr[j] * scaleF;
|
|
}
|
|
fc->Set_WgSolverExcitation_H(hvtr);
|
|
}
|
|
fclose(foo);
|
|
}
|
|
|
|
void portMesh::makeCoordSystem()
|
|
{
|
|
vtr zAxis = findFaceNormal(*(fcArray[0]));
|
|
zAxis = zAxis * (-1.0); //point inside tetr(0)
|
|
// zAxis = zAxis *(1.0); // forSMT
|
|
// forthe SMT example somehow the zAxis was (0,0,1) but from the geometry
|
|
// it should be (0,0,-1). So I flip the sign to make the right direction
|
|
coord.setz_axis(zAxis); // I always assume the zAxis is the portDirection
|
|
|
|
node* vertexNode = 0;
|
|
node* edgeNode = 0;
|
|
findNodes(&vertexNode, &edgeNode);
|
|
if(vertexNode == 0){
|
|
cerr << "ERROR: PortMesh::makeCoordSystem: " << "Could not find a vertex node forport " << portName << endl;
|
|
exit(1);
|
|
}
|
|
if(edgeNode == 0){
|
|
cerr << "ERROR: PortMesh::makeCoordSystem: " << "Could not find an edge node forport " << portName << endl;
|
|
exit(1);
|
|
}
|
|
vtr O = vertexNode->getCoord();
|
|
coord.setO(O);
|
|
|
|
vtr xAxis = edgeNode->getCoord() - O;
|
|
coord.setx_axis(xAxis);
|
|
|
|
vtr yAxis = zAxis * xAxis;
|
|
coord.sety_axis(yAxis);
|
|
vtr PortDirection = coord.getZAxis();
|
|
|
|
cout << portName << endl;
|
|
cout << "( " << PortDirection.getx() << ", " << PortDirection.gety() << ", " << PortDirection.getz() << " )" << endl;
|
|
|
|
fp_t Port_Direction_mag = sqrt(PortDirection.getx() * PortDirection.getx() + PortDirection.gety() * PortDirection.gety() + PortDirection.getz() * PortDirection.getz());
|
|
|
|
fp_t px = PortDirection.getx()/Port_Direction_mag;
|
|
fp_t py = PortDirection.gety()/Port_Direction_mag;
|
|
fp_t pz = PortDirection.getz()/Port_Direction_mag;
|
|
PortDirection_vtr.setvtr(px,py,pz);
|
|
|
|
for(int i = 0; i < faceCNT; i++)
|
|
fcArray[i]->getbcPtr()->set_PortDirection( PortDirection_vtr.getx(), PortDirection_vtr.gety(), PortDirection_vtr.getz());
|
|
}
|
|
|
|
vtr portMesh::findFaceNormal(face& Face){
|
|
int i;
|
|
|
|
vtr vv0 = Face.getNode(0)->getCoord();
|
|
vtr vv1 = Face.getNode(1)->getCoord();
|
|
vtr vv2 = Face.getNode(2)->getCoord();
|
|
vtr nvtr = (vv1 - vv0) * (vv2 - vv0);
|
|
nvtr = nvtr * 0.5;
|
|
|
|
tetra& Tetra = Face.getTetra(0);
|
|
if(Tetra == 0){
|
|
cerr << "ERROR: DgGroup::findFaceNormal: " << "Port face does not have a correct tetra link." << endl;
|
|
exit(1);
|
|
}else{
|
|
vv0.setvtr(0., 0., 0.);
|
|
for(i = 0; i < 4; i++){
|
|
vv0 = vv0 + Tetra.getNode(i)->getCoord();
|
|
}
|
|
vv0 = vv0 * 0.25;
|
|
vv0 = vv0 - vv1;
|
|
if(dotP(vv0, nvtr) > 0.)
|
|
nvtr = nvtr * (-1.);
|
|
}
|
|
nvtr = nvtr * (1.0 / nvtr.magnitude());
|
|
return nvtr;
|
|
}
|
|
|
|
void portMesh::findNodes(node** vNode, node** eNode){
|
|
int i, j;
|
|
vtr vv1;
|
|
vtr vv2;
|
|
|
|
for(i = 0; i < faceCNT; i++){
|
|
bool foundVertex = false;
|
|
bool foundEdge = false;
|
|
node* vertexNode = 0;
|
|
node* edgeNode = 0;
|
|
for(j = 0; j < 3; j++){
|
|
node* Node = fcArray[i]->getNodePtr(j);
|
|
if(Node->getPType() == VERTEX_NODE){
|
|
if(foundVertex == true){
|
|
if(foundEdge == true){
|
|
vv1 = edgeNode->getCoord() - vertexNode->getCoord();
|
|
vv2 = Node->getCoord() - vertexNode->getCoord();
|
|
if(vv2.magnitude() < vv1.magnitude()){
|
|
edgeNode = Node;
|
|
}
|
|
}else{
|
|
foundEdge = true;
|
|
edgeNode = Node;
|
|
}
|
|
}else{
|
|
vertexNode = Node;
|
|
foundVertex = true;
|
|
}
|
|
}
|
|
if(Node->getPType() == EDGE_NODE){
|
|
foundEdge = 1;
|
|
edgeNode = Node;
|
|
}
|
|
}
|
|
if((foundVertex) && (foundEdge)){
|
|
*vNode = vertexNode;
|
|
*eNode = edgeNode;
|
|
return;
|
|
}
|
|
}
|
|
}
|
|
|
|
void portMesh::writeVtk()
|
|
{
|
|
int i, j;
|
|
char hFileName[STRLEN];
|
|
char eFileName[STRLEN];
|
|
FILE* fdE;
|
|
FILE* fdH;
|
|
char cc;
|
|
int dummy;
|
|
dbl_vtr dummyVtr;
|
|
complex<double> gamma_;
|
|
double modeIMP_;
|
|
|
|
sprintf(eFileName, "%s.1E1", portName);
|
|
if((fdE = fopen(eFileName, "rb")) == 0)
|
|
{
|
|
cerr << "ERROR: PortMesh::makeRhs: " << "Could not open file " << eFileName << " forread." << endl;
|
|
exit(1);
|
|
}
|
|
|
|
fread(&dummy, sizeof(int), 1, fdE);
|
|
fread(&cc, sizeof(char), 1, fdE);
|
|
fread(&cc, sizeof(char), 1, fdE);
|
|
fread(&dummy, sizeof(int), 1, fdE);
|
|
fread(&gamma_, sizeof(complex<double>), 1, fdE);
|
|
for(i = 0; i < 3; i++)
|
|
fread(&dummyVtr, sizeof(dbl_vtr), 1, fdE);
|
|
fread(&modeIMP_, sizeof(double), 1, fdE);
|
|
|
|
sprintf(hFileName, "%s.1H1", portName);
|
|
if((fdH = fopen(hFileName, "rb")) == 0)
|
|
{
|
|
cerr << "ERROR: PortMesh::makeRhs: " << "Could not open file " << hFileName << " forread." << endl;
|
|
exit(1);
|
|
}
|
|
|
|
fread(&dummy, sizeof(int), 1, fdH);
|
|
fread(&cc, sizeof(char), 1, fdH);
|
|
fread(&cc, sizeof(char), 1, fdH);
|
|
fread(&dummy, sizeof(int), 1, fdH);
|
|
fread(&gamma_, sizeof(complex<double>), 1, fdH);
|
|
for(i = 0; i < 3; i++)
|
|
fread(&dummyVtr, sizeof(dbl_vtr), 1, fdH);
|
|
fread(&modeIMP_, sizeof(double), 1, fdH);
|
|
|
|
// fill the port field with averaged values
|
|
vtr* eField = new vtr[nodeCNT];
|
|
vtr* hField = new vtr[nodeCNT];
|
|
int* count = new int[nodeCNT];
|
|
memset(count, 0, nodeCNT*sizeof(int));
|
|
|
|
int index;
|
|
dbl_cvtr evtr;
|
|
dbl_cvtr hvtr;
|
|
vtr tmpE; // change here
|
|
vtr tmpH;
|
|
for(i = 0; i < faceCNT; i++){
|
|
for(j = 0; j < SEC_ORDER_FACE_NODES; j++){
|
|
fread(&evtr, sizeof(dbl_cvtr), 1, fdE);
|
|
fread(&hvtr, sizeof(dbl_cvtr), 1, fdH);
|
|
|
|
if(j < 3)
|
|
{
|
|
index = globToLocMap_->find(fcArray[i]->getNode(j)->getid())->second;
|
|
|
|
tmpE = coord.x_axis * (fp_t)(evtr.x.real()) + coord.y_axis * (fp_t)(evtr.y.real()) + coord.z_axis * (fp_t)(evtr.z.real());
|
|
tmpH = coord.x_axis * (fp_t)(hvtr.x.real()) + coord.y_axis * (fp_t)(hvtr.y.real()) + coord.z_axis * (fp_t)(hvtr.z.real());
|
|
|
|
eField[index] = eField[index] + tmpE;
|
|
hField[index] = hField[index] + tmpH;
|
|
(count[index])++;
|
|
}
|
|
}
|
|
}
|
|
|
|
fclose(fdE);
|
|
fclose(fdH);
|
|
|
|
for(i = 0; i < nodeCNT; i++){
|
|
eField[i] = eField[i] / static_cast<fp_t>(count[i]);
|
|
hField[i] = hField[i] / static_cast<fp_t>(count[i]);
|
|
}
|
|
|
|
node** locNodeArray = new node*[nodeCNT];
|
|
for(i = 0; i < nodeCNT; i++){
|
|
node& Node = *(ndArray[i]);
|
|
int index = globToLocMap_->find(Node.getid())->second;
|
|
locNodeArray[index] = new node(index, Node.getPType(), Node.getSingOrder(), Node.getCoord().getx(), Node.getCoord().gety(), Node.getCoord().getz());
|
|
}
|
|
|
|
face** locFaceArray = new face*[faceCNT];
|
|
for(i = 0; i < faceCNT; i++){
|
|
face& Face = *(fcArray[i]);
|
|
locFaceArray[i] = new face(Face);
|
|
locFaceArray[i]->setFace(
|
|
locNodeArray[globToLocMap_->find(Face.getNode(0)->getid())->second],
|
|
locNodeArray[globToLocMap_->find(Face.getNode(1)->getid())->second],
|
|
locNodeArray[globToLocMap_->find(Face.getNode(2)->getid())->second]);
|
|
}
|
|
|
|
VtkWriter vtkWriter(1.);
|
|
vtkWriter.writeTriUg(portName, nodeCNT, locNodeArray, faceCNT, locFaceArray, eField, hField, 1);
|
|
|
|
for(i = 0; i < nodeCNT; i++)
|
|
delete locNodeArray[i];
|
|
delete [] locNodeArray;
|
|
for(i = 0; i < faceCNT; i++)
|
|
delete locFaceArray[i];
|
|
delete [] locFaceArray;
|
|
delete [] eField;
|
|
delete [] hField;
|
|
delete [] count;
|
|
}
|
|
|
|
|
|
|
|
|
|
// -----------------------------------------------------------------------------
|
|
// portMesh::makeRHS_TEM (HFSS-style symmetry port)
|
|
//
|
|
// Builds a uniform TEM-like snapshot on the port cross-section:
|
|
//
|
|
// • Symmetry planes are applied only on the port boundary (not physical walls).
|
|
// • E is set to point across the PEC symmetry plane (normal to it, in-plane).
|
|
// • H = (1/η) (ez × E), with ez the port normal (power +S).
|
|
// • |E| is fixed to 1 V/m (no 1-W renormalization).
|
|
//
|
|
// Inputs:
|
|
// freq_Hz, eps_r : frequency [Hz] and relative permittivity (uniform medium)
|
|
// Ex_dir_* : a vector pointing normal to the PEC symmetry plane,
|
|
// lying in the port plane (if not, it will be projected).
|
|
// PortDir_* : the port normal (+S direction).
|
|
//
|
|
// Side effects / state:
|
|
// • For every face, provides 6 samples (2nd-order tri) via Set_WgSolverExcitation_E/H.
|
|
// • Stores wave (field) impedance η = sqrt(μ/ε) in modeIMP and modeImpedance.
|
|
// • Stores gamma = k = ω sqrt(μ ε) (lossless).
|
|
// • Logs the resulting power for |E|=1 V/m, i.e. P = 0.5 * (1/η) * A.
|
|
// -----------------------------------------------------------------------------
|
|
void portMesh::makeRHS_TEM(double freq_Hz, double eps_r,
|
|
double Ex_dir_x, double Ex_dir_y, double Ex_dir_z, // E direction (normal to PEC plane, in-plane)
|
|
double PortDir_x, double PortDir_y, double PortDir_z)// propagation (+S)
|
|
{
|
|
using std::complex;
|
|
|
|
std::cout << " [TEM|Sym] Fill for " << portName
|
|
<< " @ f=" << freq_Hz/1e9 << " GHz, eps_r=" << eps_r
|
|
<< " (|E| = 1 V/m)\n";
|
|
|
|
// --- 1) Physical constants and medium parameters (SI) ----------------------
|
|
const double c0 = 299792458.0;
|
|
const double mu0 = 4e-7 * M_PI;
|
|
const double eps0 = 1.0 / (mu0 * c0 * c0);
|
|
|
|
const double mu = mu0; // non-magnetic medium
|
|
const double eps = eps_r * eps0;
|
|
const double w = 2.0 * M_PI * freq_Hz;
|
|
const double k = w * std::sqrt(mu * eps); // β magnitude for lossless TEM
|
|
const double eta = std::sqrt(mu / eps); // wave impedance E/H (field impedance)
|
|
|
|
// Store in your legacy members
|
|
gamma = static_cast<fp_t>(k);
|
|
modeIMP = static_cast<fp_t>(eta);
|
|
modeImpedance = complex<double>(eta, 0.0);
|
|
|
|
// --- 2) Local orthonormal frame {ex, ey, ez} --------------------------------
|
|
// ez = port normal (+S direction)
|
|
vtr ez(PortDir_x, PortDir_y, PortDir_z);
|
|
ez = normalize(ez);
|
|
if (norm(ez) == 0.0) {
|
|
std::cerr << " [TEM|Sym] ERROR: PortDir is zero.\n";
|
|
return;
|
|
}
|
|
|
|
// ex = desired E direction (normal to PEC symmetry plane), projected to port plane
|
|
vtr ex(Ex_dir_x, Ex_dir_y, Ex_dir_z);
|
|
if (norm(ex) == 0.0) {
|
|
std::cerr << " [TEM|Sym] ERROR: E direction is zero.\n";
|
|
return;
|
|
}
|
|
// Project onto plane ⟂ ez and normalize
|
|
ex = ex - ez * dot(ex, ez);
|
|
ex = normalize(ex);
|
|
if (norm(ex) == 0.0)
|
|
{
|
|
// Given direction was (nearly) parallel to ez; pick a stable in-plane axis
|
|
std::cerr << " [TEM|Sym] WARNING: E direction ~ parallel to PortDir; choosing in-plane axis.\n";
|
|
vtr up = (std::fabs(ez.getz()) > 0.9) ? vtr(0,1,0) : vtr(0,0,1);
|
|
ex = normalize( cross(up, ez) );
|
|
if (norm(ex) == 0.0) { up = vtr(0,1,0); ex = normalize( cross(up, ez) ); }
|
|
if (norm(ex) == 0.0) {
|
|
std::cerr << " [TEM|Sym] ERROR: cannot establish in-plane E direction.\n";
|
|
return;
|
|
}
|
|
}
|
|
|
|
// ey completes right-handed frame
|
|
vtr ey = cross(ez, ex);
|
|
|
|
// --- 3) Compute port area A (sum of triangle areas) -------------------------
|
|
// OPTION A: use node* (NOT const node*) to call non-const node::getCoord()
|
|
double A_port = 0.0;
|
|
for (int f = 0; f < faceCNT; ++f) {
|
|
node* n0 = fcArray[f]->getNode(0);
|
|
node* n1 = fcArray[f]->getNode(1);
|
|
node* n2 = fcArray[f]->getNode(2);
|
|
|
|
const vtr r0(n0->getCoord().getx(), n0->getCoord().gety(), n0->getCoord().getz());
|
|
const vtr r1(n1->getCoord().getx(), n1->getCoord().gety(), n1->getCoord().getz());
|
|
const vtr r2(n2->getCoord().getx(), n2->getCoord().gety(), n2->getCoord().getz());
|
|
|
|
A_port += tri_area(r0, r1, r2);
|
|
}
|
|
if (!(A_port > 0.0)) {
|
|
std::cerr << " [TEM|Sym] ERROR: port area is zero.\n";
|
|
return;
|
|
}
|
|
|
|
// --- 4) Fix |E| = 1 V/m; H consistent with η; compute resulting power -------
|
|
const double E0 = 1.0; // requested fixed magnitude
|
|
const vtr Evec = ex * static_cast<fp_t>(E0);
|
|
const vtr Hvec = ey * static_cast<fp_t>(E0 / eta);
|
|
|
|
// For info: total time-average power P = 0.5 * (E×H)·ez * A
|
|
// Here (E×H)·ez = (E0^2)/η and is uniform over the port
|
|
const double P_result = 0.5 * ( (E0 * E0) / eta ) * A_port;
|
|
|
|
// --- 5) Provide 6 samples per face (uniform, so same vector at each sample) -
|
|
vtr E6[6], H6[6];
|
|
for (int f = 0; f < faceCNT; ++f)
|
|
{
|
|
E6[0] = Evec; H6[0] = Hvec; // vertex 0
|
|
E6[1] = Evec; H6[1] = Hvec; // vertex 1
|
|
E6[2] = Evec; H6[2] = Hvec; // vertex 2
|
|
E6[3] = Evec; H6[3] = Hvec; // midside (0,1)
|
|
E6[4] = Evec; H6[4] = Hvec; // midside (1,2)
|
|
E6[5] = Evec; H6[5] = Hvec; // midside (2,0)
|
|
|
|
fcArray[f]->Set_WgSolverExcitation_E(E6);
|
|
fcArray[f]->Set_WgSolverExcitation_H(H6);
|
|
|
|
cout << " Face " << f << ": E = ("
|
|
<< Evec.getx() << ", " << Evec.gety() << ", " << Evec.getz() << ") V/m, |E|=" << norm(Evec) << " V/m\n";
|
|
}
|
|
|
|
// --- 6) Log summary ---------------------------------------------------------
|
|
std::cout << " [TEM|Sym] η = " << eta
|
|
<< " Ω (E/H), A = " << A_port << " m^2, |E| = " << E0 << " V/m"
|
|
<< ", resulting P = " << P_result << " W\n";
|
|
} |