#include <stdio.h>
#include <stdlib.h>
#include "libsx.h"
#include "globals.h"
#include "grid.h"
#include "graph.h"
#include "kraftwerk.h"
#include "post.h"
#include "module.h"
#include "geodata.h"
#define CANT_FIND 0
#define FOUND_IT 1
typedef struct {
int seg;
int n;
int *nodes;
} segmentsOnNodes;
segmentsOnNodes *segnod;
int nseg;
float xmin_orig;
float xmax_orig;
float ymin_orig;
float ymax_orig;
int positionFilePointer(FILE *stream, char *tag);
void readEntry(char *file, char *tag, void (*readingFunction)() );
void addSolutionMenus();
int thisType;
int eqDataItem;
void readMesh(FILE *stream)
{
int i;
if(stream != NULL){
currentSolution = 0;
SetToggleState(solButt[currentSolution], TRUE);
NoOfMeshes = 1;
fscanf(stream, "%i%i", &mySol[currentSolution].nnode, &mySol[currentSolution].elnum);
if(mySol[currentSolution].elnum == 0)
errorHandler("Invalid data file!");
if(allocateMemory(currentSolution, mySol[currentSolution].elnum,
mySol[currentSolution].nnode) == 0){
for( i = 0; i < mySol[currentSolution].nnode; i++)
fscanf(stream,"%f%f", &mySol[currentSolution].coord[i][0] ,&mySol[currentSolution].coord[i][1]);
for( i = 0; i < mySol[currentSolution].elnum; i++)
fscanf(stream,"%i%i%i", &mySol[currentSolution].nodpek[i][0],&mySol[currentSolution].nodpek[i][1],&mySol[currentSolution].nodpek[i][2]) ;
/*
* Renumber node numbers so they start from zero
*
*/
for(i = 0; i < mySol[currentSolution].elnum; i++)
{
mySol[currentSolution].nodpek[i][0] -= 1;
mySol[currentSolution].nodpek[i][1] -= 1;
mySol[currentSolution].nodpek[i][2] -= 1;
}
}
}
else
errorHandler("Unable to open file");
}
void new_readPoints(FILE *stream)
{
int i, j, n;
float f1, f2;
PTS_NODE *p;
fscanf(stream, "%i", &n);
for(i = 0; i < n; i++){
fscanf(stream, "%i %f %f", &j, &f1, &f2);
p = add_pnt(f1, f2);
p->id = j;
}
}
void new_readSegments(FILE *stream)
{
PTS_NODE *pts1, *pts2;
GEO_NODE *b;
BND_NODE *s;
int n, i, id;
int p1, p2;
fscanf(stream, "%i ", &n);
for(i = 0; i < n; i++){
fscanf(stream, "%i %i %i ", &id, &p1, &p2);
/* Now find the pointers associated with p1,p2 */
pts1 = get_pts_pointer(p1);
pts2 = get_pts_pointer(p2);
FOR_ALL_BOUNDARIES(b){
FOR_ALL_SEGMENTS(b, s){
if(s->seg->id == id){
s->seg->pts1 = pts1;
s->seg->pts2 = pts2;
break;
}
}
}
}
}
void new_readBoundary(FILE *stream)
{
int j, id, n, dummy;
GEO_NODE *g;
fscanf(stream, "%i", &dummy);
fscanf(stream, "%i %i", &j, &n);
g = new_bnd();
g->type = thisType;
for(j = 0; j < n; j++){
fscanf(stream, "%i ", &id);
add_seg2(id, g);
}
}
void readDomain(FILE *stream)
{
fscanf(stream,"%f %f %f %f", &xmin, &xmax, &ymin, &ymax);
xmin_orig = xmin;
xmax_orig = xmax;
ymin_orig = ymin;
ymax_orig = ymax;
}
void readState(FILE *stream)
{
fscanf(stream,"%i", &state);
}
void readModuleID(FILE *stream)
{
int id;
fscanf(stream,"%i", &id);
module[id].state = state;
setModule(id);
}
void readSolution(FILE *stream)
{
int el, imesh, i;
int slask;
currentSolution = 0;
SetToggleState(solButt[currentSolution], TRUE);
/* free old pointers */
for( imesh = 0; imesh < NoOfMeshes; imesh++){
free(mySol[imesh].nodpek);
free(mySol[imesh].coord);
free(mySol[imesh].u);
free(mySol[imesh].f);
free(mySol[imesh].vecX);
free(mySol[imesh].vecY);
}
if(stream != NULL){
fscanf(stream, "%i", &NoOfMeshes);
for(imesh = 0; imesh < NoOfMeshes; imesh++){
fscanf(stream, "%i", &mySol[imesh].elnum);
fscanf(stream, "%i", &mySol[imesh].nnode);
if(allocateMemory(imesh, mySol[imesh].elnum, mySol[imesh].nnode) == 0){
for(el = 0; el < mySol[imesh].elnum; el++)
fscanf(stream, "%i%i%i", &mySol[imesh].nodpek[el][0], &mySol[imesh].nodpek[el][1],
&mySol[imesh].nodpek[el][2]);
for(el = 0; el < mySol[imesh].elnum; el++)
fscanf(stream, "%i%i%i", &slask, &slask, &slask);
for(i = 0; i < mySol[imesh].nnode; i++){
fscanf(stream, "%f", &mySol[imesh].coord[i][0]);
fscanf(stream, "%f", &mySol[imesh].coord[i][1]);
}
for(i = 0; i < mySol[imesh].nnode; i++)
fscanf(stream, "%f", &mySol[imesh].u[i]);
for(i = 0; i < mySol[imesh].nnode; i++)
fscanf(stream, "%f%f", &mySol[imesh].vecX[i], &mySol[imesh].vecY[i]);
for(el = 0; el < mySol[imesh].elnum; el++){
mySol[imesh].nodpek[el][0] -= 1;
mySol[imesh].nodpek[el][1] -= 1;
mySol[imesh].nodpek[el][2] -= 1;
}
} /* Next mesh */
}
}
else
errorHandler("Unable to open the solution file");
}
void readSolution2(FILE *stream)
{
int el, imesh, i, j;
int slask, ioff, ntot;
char dummy;
currentSolution = 0;
SetToggleState(solButt[currentSolution], TRUE);
/* free old pointers */
for( imesh = 0; imesh < NoOfMeshes; imesh++){
free(mySol[imesh].nodpek);
free(mySol[imesh].coord);
free(mySol[imesh].u);
free(mySol[imesh].f);
free(mySol[imesh].vecX);
free(mySol[imesh].vecY);
}
if(stream != NULL){
fscanf(stream, "%i", &NoOfMeshes);
for(imesh = 0; imesh < NoOfMeshes; imesh++){
fscanf(stream, "%i", &mySol[imesh].elnum);
fscanf(stream, "%i", &mySol[imesh].nnode);
fscanf(stream, "%i %i ", &mySol[imesh].n, &ntot);
if(allocateMemory(imesh, mySol[imesh].elnum, mySol[imesh].nnode) == 0){
for(el = 0; el < mySol[imesh].elnum; el++)
fscanf(stream, "%i%i%i", &mySol[imesh].nodpek[el][0], &mySol[imesh].nodpek[el][1],
&mySol[imesh].nodpek[el][2]);
for(el = 0; el < mySol[imesh].elnum; el++)
fscanf(stream, "%i%i%i", &slask, &slask, &slask);
for(i = 0; i < mySol[imesh].nnode; i++){
fscanf(stream, "%f ", &mySol[imesh].coord[i][0]);
fscanf(stream, "%f ", &mySol[imesh].coord[i][1]);
}
mySol[imesh].f = malloc(ntot * mySol[imesh].nnode * sizeof(float));
ioff = 0;
for(j = 0; j < mySol[imesh].n; j++){
mySol[imesh].ptr[j] = ioff;
fscanf(stream, "%[^;]", mySol[imesh].name[j]);
fscanf(stream, "%c", &dummy);
fscanf(stream, " %i ", &mySol[imesh].type[j]);
printf(" %s\n", mySol[imesh].name[j]);
for(i = 0; i < mySol[imesh].nnode; i++){
if(mySol[imesh].type[j] == 1){
fscanf(stream, "%f ", &mySol[imesh].f[i + ioff]);
}
else{
fscanf(stream, "%f %f ", &mySol[imesh].f[i + ioff],
&mySol[imesh].f[i + ioff + mySol[imesh].nnode]);
}
}
if(mySol[imesh].type[j] == 1)
ioff += mySol[imesh].nnode;
else
ioff += 2*mySol[imesh].nnode;
}
for(el = 0; el < mySol[imesh].elnum; el++){
mySol[imesh].nodpek[el][0] -= 1;
mySol[imesh].nodpek[el][1] -= 1;
mySol[imesh].nodpek[el][2] -= 1;
}
} /* Next mesh */
}
}
else
errorHandler("Unable to open the solution file");
}
void readEqData(FILE *stream)
{
int i, j;
float x, y;
int item;
char dummy;
item = eqDataItem;
fscanf(stream, "%i ", &module[currentModule].data[item].n);
fscanf(stream, "%i ", &module[currentModule].dataGroup[item].type);
if(module[currentModule].dataGroup[item].type > 0){
for(j = 0; j < module[currentModule].data[item].n; j++){
if(module[currentModule].dataGroup[item].type == 2) /* elementwise */
fscanf(stream, "%i ", &module[currentModule].data[item].value[j].el);
else{ /* pointwise */
fscanf(stream, "%f %f ", &x, &y);
phys2pix(&module[currentModule].data[item].value[j].x,
&module[currentModule].data[item].value[j].y, x, y);
}
for(i = 0; i < module[currentModule].dataGroup[item].n; i++){
fscanf(stream, "%f ", &module[currentModule].data[item].value[j].f[i]);
}
}
}
else{ /* analytical */
for(i = 0; i < module[currentModule].dataGroup[item].n; i++){
fscanf(stream, "%[^;]", module[currentModule].data[item].expr[i].text);
fscanf(stream, "%c", &dummy);
}
}
}
void readNodesOnSegments(FILE *stream)
{
int i, j, n, node;
fscanf(stream, "%i", &nseg);
free(segnod);
segnod = malloc( nseg * sizeof(segmentsOnNodes));
for(i = 0; i < nseg; i++){
fscanf(stream, "%i", &n);
segnod[i].seg = n;
fscanf(stream, "%i", &n);
segnod[i].n = n;
segnod[i].nodes = malloc(n * sizeof(int));
for(j = 0; j < n; j++){
fscanf(stream, "%i", &node);
segnod[i].nodes[j] = node;
}
}
}
void readBCData(FILE *stream)
{
int ibc;
int k;
for(nvar = 0; nvar < module[currentModule].nvar; nvar++){
fscanf(stream, "%i ", &module[currentModule].bc[nvar].n);
for(ibc = 0; ibc < module[currentModule].bc[nvar].n; ibc++){
fscanf(stream, "%i %c ", &module[currentModule].bc[nvar].segment[ibc],
&module[currentModule].bc[nvar].type[ibc]);
fscanf(stream, "%f ", &module[currentModule].bc[nvar].value[ibc]);
for(k = 0; k < module[currentModule].nbc; k++){
if(module[currentModule].bc[nvar].type[ibc] ==
module[currentModule].bc[k].tag[nvar]){
module[currentModule].bc[nvar].id[ibc] = k;
}
}
}
}
}
void readNvar(FILE *stream)
{
fscanf(stream, "%i", &nvar);
}
void loadFromFile(char *file)
{
int i, j, n, x, y;
int active[20], oldstate;
float xp, yp;
PTS_NODE *p;
oldstate = state;
nvar = 0;
state = HAVE_NOTHING;
clear_geometry();
readEntry(file, "DOMAIN", readDomain);
readEntry(file, "STATE", readState);
readEntry(file, "MODULE ID", readModuleID);
readEntry(file, "POINTS", new_readPoints);
readEntry(file, "MESH", readMesh);
readEntry(file, "NODES ON SEGMENTS", readNodesOnSegments);
thisType = INNER;
readEntry(file, "INNER BOUNDARY", new_readBoundary);
thisType = OUTER;
readEntry(file, "OUTER BOUNDARY", new_readBoundary);
thisType = INTERNAL;
readEntry(file, "INTERNAL BOUNDARY", new_readBoundary);
readEntry(file, "SEGMENTS", new_readSegments);
nvar = 0;
readEntry(file, "BOUNDARY CONDITIONS", readBCData);
for(i = 0; i < 20; i++)
active[i] = FALSE;
if(oldstate == HAVE_SOLUTION){
for(i = 0; i < mySol[currentSolution].n + 1; i++)
active[i] = GetMenuItemChecked(plotMenu.item[i+5]);
}
else
active[0] = TRUE;
rmSolutionMenus();
readEntry(file, "SOLUTION", readSolution2);
addSolutionMenus();
for(i = 0; i < mySol[currentSolution].n + 1; i++)
SetMenuItemChecked(plotMenu.item[i+5], active[i]);
getDrawingScale();
for( eqDataItem = 0; eqDataItem < module[currentModule].ngroups; eqDataItem++)
readEntry(file, module[currentModule].dataGroup[eqDataItem].tag, readEqData);
FOR_ALL_POINTS(p){
xp = p->physicalX;
yp = p->physicalY;
phys2pix(&x, &y, xp, yp);
p->x = x;
p->y = y;
}
}
/*
* Below are functions needed to save things to file
*
*/
void writeNodesOnSegments(FILE *stream)
{
int i, j;
if( nseg > 0){
fprintf(stream, "%s\n", "NODES ON SEGMENTS");
fprintf(stream, "%i\n", nseg);
for(i = 0; i < nseg; i++){
fprintf(stream, "%i %i\n", segnod[i].seg, segnod[i].n);
for(j = 0; j < segnod[i].n; j++){
fprintf(stream, "%i\n", segnod[i].nodes[j]);
}
}
}
}
void new_writeMeshData(FILE *stream)
{
/*
* Dumps the boundary to the scratch file needed for meshgeneration
*
*/
int i, j, seg, ibc;
int el, imesh;
float x, y;
BND_NODE *s;
GEO_NODE *g;
PTS_NODE *p;
/* Convert the X11-integer coordinates to physical coordinates */
if(state <= HAVE_MESHSIZE)
getDrawingScale();
/* Put to file */
if(stream != NULL){
/* Put a number on each point (needed for the fortran prog) */
/*
printf("RENUMBERING ALL GEOMETRY STUFF ***************\n");
i = 1;
FOR_ALL_POINTS(p){
p->id = i;
i++;
}
i = 1;
FOR_ALL_BOUNDARIES( g ){
FOR_ALL_SEGMENTS(g, s){
s->seg->id = i;
i++;
}
}
*/
fprintf(stream, "%s\n","POINTS");
fprintf(stream, "%i\n", number_of_points());
FOR_ALL_POINTS(p)
fprintf(stream, "%i %f %f\n", p->id , p->physicalX, p->physicalY);
i = 1;
FOR_ALL_BOUNDARIES( g ){
fprintf(stream,"%s\n","SEGMENTS");
fprintf(stream,"%i\n", number_of_segments(g));
FOR_ALL_SEGMENTS(g, s)
fprintf(stream,"%i %i %i \n", s->seg->id, s->seg->pts1->id, s->seg->pts2->id);
if(g->type == INNER)
fprintf(stream,"%s\n","INNER BOUNDARY");
else if(g->type == OUTER)
fprintf(stream,"%s\n","OUTER BOUNDARY");
else
fprintf(stream,"%s\n","INTERNAL BOUNDARY");
fprintf(stream,"%i\n", 1);
fprintf(stream,"%i %i \n", i, number_of_segments(g));
FOR_ALL_SEGMENTS(g, s)
fprintf(stream, "%i \n", s->seg->id);
i++;
}
if(ih > 0){
fprintf(stream,"%s\n", "SIZE");
fprintf(stream,"%i\n", ih);
for(i = 0; i < ih; i++){
pix2phys(hx[i], hy[i], &x, &y);
fprintf(stream,"%f %f %f\n", x, y, h[i]);
}
}
if(state >= HAVE_MESH){
imesh = 0;
fprintf(stream,"%s\n", "MESH DATA");
fprintf(stream, "%i %i\n", mySol[imesh].nnode, mySol[imesh].elnum);
for(i = 0; i < mySol[imesh].nnode; i++){
fprintf(stream, "%f\n", mySol[imesh].coord[i][0]);
fprintf(stream, "%f\n", mySol[imesh].coord[i][1]);
}
for(el = 0; el < mySol[imesh].elnum; el++)
fprintf(stream, "%i %i %i\n", 1 + mySol[imesh].nodpek[el][0],
1 + mySol[imesh].nodpek[el][1],
1 + mySol[imesh].nodpek[el][2]);
writeNodesOnSegments(stream);
}
}
else
errorHandler("Unable to create area.scratch");
}
void writeSolverData(FILE *stream)
{
fprintf(stream, "%s\n", "REFINEMENT DATA");
fprintf(stream, "%f %i '%s'\n", solverTol, solverLevels, solverNorm);
}
void writeEqData(FILE *stream)
{
int item, i, j, ibc;
float x, y;
for(item = 0; item < module[currentModule].ngroups; item++){
fprintf(stream, "%s\n", module[currentModule].dataGroup[item].tag);
fprintf(stream, "%i\n", module[currentModule].data[item].n);
fprintf(stream, "%i\n", module[currentModule].dataGroup[item].type);
if(module[currentModule].dataGroup[item].type == 0){
for(i = 0; i < module[currentModule].dataGroup[item].n; i++)
fprintf(stream, "%s;", module[currentModule].data[item].expr[i].text);
fprintf(stream, "\n");
}
else{
for(j = 0; j < module[currentModule].data[item].n; j++){
if(module[currentModule].dataGroup[item].type == 1){
pix2phys(module[currentModule].data[item].value[j].x,
module[currentModule].data[item].value[j].y , &x, &y);
fprintf(stream, "%f %f ", x, y);
}
else{
fprintf(stream, "%i ", module[currentModule].data[item].value[j].el);
}
for(i = 0; i < module[currentModule].dataGroup[item].n; i++)
fprintf(stream, "%f ", module[currentModule].data[item].value[j].f[i]);
fprintf(stream, "\n");
}
}
}
/*
for( item = 0; item < numberOfEqData; item++){
fprintf(stream, "%s\n", cddata[item].label);
fprintf(stream, "%i\n", cddata[item].n);
for(j = 0; j < cddata[item].n; j++){
pix2phys(cddata[item].value[j].x, cddata[item].value[j].y, &x, &y);
fprintf(stream, "%f %f ", x, y);
for(i = 0; i < cddata[item].nvalues; i++){
fprintf(stream, "%f ", cddata[item].value[j].f[i]);
}
fprintf(stream, "%s\n"," ");
}
}
*/
fprintf(stream, "%s\n", "NUMBER OF VARIABLES");
fprintf(stream, "%i\n", module[currentModule].nvar);
fprintf(stream, "%s\n", "BOUNDARY CONDITIONS");
for(j = 0; j < module[currentModule].nvar; j++){
fprintf(stream, "%i\n", module[currentModule].bc[j].n);
for(ibc = 0; ibc < module[currentModule].bc[j].n; ibc++){
fprintf(stream, "%i %c\n", module[currentModule].bc[j].segment[ibc],
module[currentModule].bc[j].type[ibc]);
fprintf(stream, "%f\n", module[currentModule].bc[j].value[ibc]);
}
}
}
void writeSolution(FILE *stream)
{
int el, imesh, i, j, ioff, ntot;
ntot = 0;
for(i = 0; i < NoOfMeshes; i++){
for(j = 0; j < mySol[i].n; j++)
ntot += mySol[i].type[j];
}
if(stream != NULL){
fprintf(stream, "%s\n", "SOLUTION");
fprintf(stream, "%i\n", NoOfMeshes);
for(imesh = 0; imesh < NoOfMeshes; imesh++){
fprintf(stream, "%i\n", mySol[imesh].elnum);
fprintf(stream, "%i\n", mySol[imesh].nnode);
fprintf(stream, "%i %i\n ", mySol[imesh].n, ntot);
for(el = 0; el < mySol[imesh].elnum; el++)
fprintf(stream, "%i %i %i\n", 1 + mySol[imesh].nodpek[el][0],
1 + mySol[imesh].nodpek[el][1],
1 + mySol[imesh].nodpek[el][2]);
for(el = 0; el < mySol[imesh].elnum; el++)
fprintf(stream, "%i %i %i\n", 1, 1, 1);
for(i = 0; i < mySol[imesh].nnode; i++){
fprintf(stream, "%f\n", mySol[imesh].coord[i][0]);
fprintf(stream, "%f\n", mySol[imesh].coord[i][1]);
}
ioff = 0;
for(j = 0; j < mySol[imesh].n; j++){
mySol[imesh].ptr[j] = ioff;
fprintf(stream, "%s;\n", mySol[imesh].name[j]);
fprintf(stream, " %i\n ",mySol[imesh].type[j]);
for(i = 0; i < mySol[imesh].nnode; i++){
if(mySol[imesh].type[j] == 1){
fprintf(stream, "%f\n", mySol[imesh].f[i + ioff]);
}
else{
fprintf(stream, "%f %f \n", mySol[imesh].f[i + ioff],
mySol[imesh].f[i + ioff + mySol[imesh].nnode]);
}
}
if(mySol[imesh].type[j] == 1)
ioff += mySol[imesh].nnode;
else
ioff += 2*mySol[imesh].nnode;
}
} /* Next mesh */
}
else
errorHandler("Unable to open the solution file");
}
void saveToFile(char *file)
{
FILE *stream;
getDrawingScale();
stream = fopen(file, "w");
if(stream != NULL){
fprintf(stream, "#---------------------------------------------- \n");
fprintf(stream, "# Data saved by FEMLAB\n");
fprintf(stream, "# FEMLAB is brought to you by:\n");
fprintf(stream, "# Dep. of Mathematics\n");
fprintf(stream, "# Chalmers University of Technology\n");
fprintf(stream, "# Gothenburg Sweden\n");
fprintf(stream, "# (femlab@math.chalmers.se)\n");
fprintf(stream, "# \n");
fprintf(stream, "# FEMLAB, the choice of generation X\n");
fprintf(stream, "#---------------------------------------------- \n");
fprintf(stream, "%s\n", "MODULE ID");
fprintf(stream, "%i\n", module[currentModule].id);
fprintf(stream, "%s\n", module[currentModule].name);
fprintf(stream, "%s\n", "STATE");
fprintf(stream, "%i\n", state);
fprintf(stream, "%s\n", "DOMAIN");
fprintf(stream,"%f %f %f %f\n", xmin, xmax, ymin, ymax);
if(state >= HAVE_BOUNDARY){
new_writeMeshData(stream);
}
if(state >= HAVE_EQDATA)
writeEqData(stream);
if(state >= HAVE_SOLVDATA)
writeSolverData(stream);
if(state >= HAVE_SOLUTION)
writeSolution(stream);
fprintf(stream, "%s\n","END OF DATA");
fclose(stream);
}
else
errorHandler("Unable to create file");
}
void readEntry(char *file, char *tag, void (*readingFunction)() )
{
FILE *stream;
stream = fopen(file, "r");
if(stream != NULL){
while( positionFilePointer(stream, tag) != CANT_FIND){
printf("reading: %s\n", tag);
readingFunction(stream);
}
fclose(stream);
}
}
int positionFilePointer(FILE *stream, char *tag)
{
char line[80];
fgets(line, 80, stream);
while( strstr(line, tag) == NULL && strstr(line, "END") == NULL){
fgets(line, 80, stream);
}
if(strstr(line, "END") != NULL ){
return CANT_FIND;
}
else{
/*
We have a problem with keywords SEGMENTS and NODES ON SEGMENTS
It is not possible to detect which of these tags we have found
The following quick hack tries to resolve this bug
*/
if(strstr(line,"NODES ON SEGMENTS") != NULL &&
strstr(tag, "NODES ON SEGMENTS")==NULL)
return CANT_FIND;
else
return FOUND_IT;
}
}
syntax highlighted by Code2HTML, v. 0.9.1