/*****************************************************************************1*2* Elmer, A Finite Element Software for Multiphysical Problems3*4* Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland5*6* This library is free software; you can redistribute it and/or7* modify it under the terms of the GNU Lesser General Public8* License as published by the Free Software Foundation; either9* version 2.1 of the License, or (at your option) any later version.10*11* This library is distributed in the hope that it will be useful,12* but WITHOUT ANY WARRANTY; without even the implied warranty of13* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU14* Lesser General Public License for more details.15*16* You should have received a copy of the GNU Lesser General Public17* License along with this library (in file ../LGPL-2.1); if not, write18* to the Free Software Foundation, Inc., 51 Franklin Street,19* Fifth Floor, Boston, MA 02110-1301 USA20*21*****************************************************************************/2223/*******************************************************************************24*25* MATC main module.26*27*******************************************************************************28*29* Author: Juha Ruokolainen30*31* Address: CSC - IT Center for Science Ltd.32* Keilaranta 14, P.O. BOX 40533* 02101 Espoo, Finland34* Tel. +358 0 457 272335* Telefax: +358 0 457 230236* EMail: [email protected]37*38* Date: 30 May 199639*40* Modified by:41*42* Date of modification:43*44******************************************************************************/45/***********************************************************************46|47| MATC - Last Edited 9. 8. 198848|49***********************************************************************/5051/*======================================================================52|Syntax of the manual pages:53|54|FUNCTION NAME(...) params ...55|56$ usage of the function and type of the parameters57? explain the effects of the function58= return value and the type of value if not of type int59@ globals effected directly by this routine60! current known bugs or limitations61& functions called by this function62~ these functions may interest you as an alternative function or63| because they control this function somehow64^=====================================================================*/656667/*68* $Id: matc.c,v 1.7 2007/06/08 08:12:17 jpr Exp $69*70* $Log: matc.c,v $71* Revision 1.7 2007/06/08 08:12:17 jpr72* *** empty log message ***73*74* Revision 1.6 2006/02/07 10:21:42 jpr75* Changed visibility of some variables to local scope.76*77* Revision 1.5 2006/02/02 06:54:44 jpr78* small formatting changes.79*80* Revision 1.3 2005/08/25 13:44:22 vierinen81* windoze stuff82*83* Revision 1.2 2005/05/27 12:26:20 vierinen84* changed header install location85*86* Revision 1.1.1.1 2005/04/14 13:29:14 vierinen87* initial matc automake package88*89* Revision 1.2 1998/08/01 12:34:48 jpr90*91* Added Id, started Log.92*93*94*/9596#define MODULE_MATC97#include "elmer/matc.h"98#include "str.h"99#include "../config.h"100101#ifdef DEBUG102static FILE *fplog;103static int tot;104#pragma omp threadprivate (fplog, tot)105#endif106/*======================================================================107? main program, initialize few constants and go for it.108^=====================================================================*/109void mtc_init( FILE *input_file, FILE *output_file, FILE *error_file )110{111VARIABLE *ptr;112113char str[256];114115int i; /* i'm getting tired with all these i's */116117static char *evalHelp =118{119"eval( str )\n\n"120"Evaluate content variable. Another form of this command is @str.\n"121};122123static char *sourceHelp =124{125"source( name )\n\n"126"Execute commands from file given name.\n"127};128129static char *helpHelp =130{131"help or help(\"symbol\")\n\n"132"First form of the command gives list of available commands.\n"133"Second form gives help on specific routine.\n"134};135136#ifdef _OPENMP137/* Allocate listheaders for each thread separately */138#pragma omp parallel139{140/* Do malloc and initialize listheaders */141listheaders = (LIST *) malloc(sizeof(LIST)*MAX_HEADERS);142/* memory allocations */143listheaders[ALLOCATIONS].next = NULL;144listheaders[ALLOCATIONS].name = "Allocations";145/* global CONSTANTS */146listheaders[CONSTANTS].next = NULL;147listheaders[CONSTANTS].name = "Constants";148/* global VARIABLES */149listheaders[VARIABLES].next = NULL;150listheaders[VARIABLES].name = "Currently defined VARIABLES";151/* internal commands */152listheaders[COMMANDS].next = NULL;153listheaders[COMMANDS].name = "Builtin Functions";154/* user defined functions */155listheaders[FUNCTIONS].next = NULL;156listheaders[FUNCTIONS].name = "User Functions";157}158#endif /* _OPENMP */159160#ifdef DEBUG161fplog = fopen("matcdbg","w");162#endif163ALLOC_HEAD = (LIST *)NULL;164165/*166* input & output & error streams167*/168math_in = input_file;169math_err = error_file;170math_out = output_file;171172mtr_com_init(); /* initialize matrix handling commands */173var_com_init(); /* "" VARIABLE "" "" */174fnc_com_init(); /* "" function handling commands */175fil_com_init(); /* "" file handling commands */176gra_com_init(); /* "" graphics commands */177str_com_init(); /* "" string handling */178179/*180* and few others.181*/182com_init( "eval" , FALSE, FALSE, com_apply, 1, 1, evalHelp );183com_init( "source" , FALSE, FALSE, com_source, 1, 1, sourceHelp );184com_init( "help" , FALSE, FALSE, com_help , 0, 1, helpHelp );185com_init( "quit" , FALSE, FALSE, com_quit , 0, 0, "quit\n" );186com_init( "exit" , FALSE, FALSE, com_quit , 0, 0, "exit\n" );187188/*189* these constants will always be there for you.190*/191ptr = const_new("true", TYPE_DOUBLE, 1, 1);192M(ptr,0,0) = 1.0;193194ptr = const_new("false", TYPE_DOUBLE, 1, 1);195M(ptr,0,0) = 0.0;196197ptr = const_new("stdin", TYPE_DOUBLE, 1, 1);198M(ptr,0,0) = 0;199200ptr = const_new("stdout", TYPE_DOUBLE, 1, 1);201M(ptr,0,0) = 1;202203ptr = const_new("stderr", TYPE_DOUBLE, 1, 1);204M(ptr,0,0) = 2;205206ptr = const_new("pi", TYPE_DOUBLE, 1, 1);207M(ptr,0,0) = 2*acos(0.0);208209#if 0210/*211* trap INTERRUPT and Floating Point Exception signals212*/213signal(SIGFPE, sig_trap);214215sprintf( str, "%s/lib/mc.ini", getenv("ELMER_POST_HOME") );216217if ( (math_in = fopen( str, "r" ) ) != (FILE *)NULL)218{219doread();220fclose( math_in );221}222223/*224* and finally standard input.225*/226math_in = stdin;227228doread();229230var_free();231com_free();232fnc_free();233const_free();234235mem_free_all();236237#ifdef DEBUG238fclose(fplog);239#endif240#endif241242return; /* done */243}244245char * mtc_domath( char *str )246{247VARIABLE *headsave; /* this should not be here */248249jmp_buf jmp, *savejmp; /* save program context */250251void (*sigfunc)() = (void (*)())signal( SIGINT, sig_trap );252253if ( !str || !*str )254{255str = (char *)doread();256signal( SIGINT, sigfunc );257return math_out_str;258}259260savejmp = jmpbuf;261jmpbuf = &jmp;262263#ifdef DEBUG264fprintf( stderr, "got [%s]\n", str );265#endif266if ( math_out_str ) math_out_str[0] = '\0';267math_out_count = 0;268269/*270* try it271*/272if (*str != '\0')273{274ALLOC_HEAD = (LIST *)NULL;275headsave = (VARIABLE *)VAR_HEAD;276277/*278* normal return takes branch 1,279* error() takes branch 2,280* quit() takes branch 3.281*/282switch (setjmp(*jmpbuf))283{284case 0:285(void)doit( str );286longjmp(*jmpbuf, 1);287break;288289case 1:290break;291292case 2:293VAR_HEAD = (LIST *)headsave;294break;295296case 3:297break;298}299}300301jmpbuf = savejmp;302303signal( SIGINT, sigfunc );304305return math_out_str;306}307308char *doread(void)309/*======================================================================310? doread() is really the main loop of this program. Function reads311| it's input as strings and gives them to function doit() for312| execution. setjmp() function is used for error recovery.313|314| Memory allocated during the lifetime of this function is315| collected to a list represented by the global VARIABLE316| ALLOCLIST *alloc_list. If function error() is called, this317| list is used to deallocate memory. Normally (well I certainly318| hope so) functions which allocate memory deallocate it themselves.319|320| Program stays in this function until an end of file -condition321| is reached or exit- or quit-commands are given.322|323@ jmp_buf *jmpbuf, ALLOC_LIST *alloc_list324& ALLOCMEM, FREEMEM, setjmp(), longjmp()325~ doit(), quit(), error()326^=====================================================================*/327{328VARIABLE *headsave; /* this should not be here */329330jmp_buf jmp, *savejmp; /* save program context */331332char *p, *q; /* buffer for input stream */333334savejmp = jmpbuf;335jmpbuf = &jmp;336337if ( math_out_str ) math_out_str[0] = '\0';338math_out_count = 0;339340p = q = ALLOCMEM(4096);341/*342* try it343*/344while(dogets(p, PMODE_MAIN))345{346if (*p != '\0')347{348ALLOC_HEAD = (LIST *)NULL;349headsave = (VARIABLE *)VAR_HEAD;350351/*352* normal return takes branch 1,353* error() takes branch 2,354* quit() takes branch 3.355*/356switch (setjmp(*jmpbuf))357{358case 0:359(void)doit(p);360longjmp(*jmpbuf, 1);361break;362363case 1:364break;365366case 2:367VAR_HEAD = (LIST *)headsave;368break;369370case 3:371goto ret;372break;373}374}375}376377ret:378379jmpbuf = savejmp;380381FREEMEM(q);382383return math_out_str;384}385386VARIABLE *com_quit(void)387/*======================================================================388? Quit current doread entry by longjumping back to it (nasty).389& longjmp390~ doread391^=====================================================================*/392{393longjmp(*jmpbuf, 3);394395return (VARIABLE *)NULL; /* won't be executed (hopefully) */396}397398int dogets(char *buff, char *prompt)399/*======================================================================400? Get line from input stream. If both input & output streams are401| connected to terminal, this function gives user one of three402| (default) prompts:403|404| MATC>405| - normal prompt (PMODE_MAIN)406| ....>407| - begin end- block is being defined (PMODE_BLOCK)408| ####> (PMODE_CONT)409| - user has given a #-sign as a last character of410| previous line, this line will be concatenated with it411|412| If current comment character is found from input stream, the413| line after this character is discarded. Likewise if current414| system command character is found, the rest of the line is415| passed to system()-call.416|417= line got -> TRUE, EOF -> FALSE418! There should be a way to get an echo when reading from file.419& fprintf(), isatty(), fileno(), strlen(), fgets(), system()420^=====================================================================*/421{422char *ptr = buff, *p; /* Can't get rid of these. */423424if ( !math_in ) return FALSE;425426/*427Try figuring out if input & output streams are428terminals, if they both are, give user a prompt.429*/430if (isatty(fileno(math_in)) && isatty(fileno(math_out)))431PrintOut( "%s", prompt );432433/*434i'm not in the mood to explain this.435*/436*ptr++ = ' ';437438/*439Go for it.440*/441while((ptr = fgets(ptr, 256, math_in)) != NULL)442{443444ptr[strlen(ptr)-1] = '\0';445446/*447* Check if the user wants to continue with this line.448*/449while(ptr[strlen(ptr)-1] == '\\')450{451ptr += strlen(ptr) - 1;452dogets(ptr, PMODE_CONT);453}454455/*456* if there is only spaces in this line,457* don't bother returning it, instead458* let's read afresh, otherwise return.459*/460p = ptr; while(isspace(*p)) p++;461462if (*p != '\0') /* GOOD EXIT HERE */463{464#if 0465/*466* Look for the system character, if found467* pass rest of the line to system()-call468*/469for(p = buff; *p; p++)470{471switch(*p)472{473case SYSTEM:474system(p + 1);475PrintOut("\n");476*p = '\0'; p--;477break;478}479}480#endif481if (*buff != '\0')482return TRUE; /* OR IF WE ARE HONEST, IT'S HERE */483}484485/*486if it's terminal give a prompt.487*/488if (isatty(fileno(math_in)) && isatty(fileno(math_out)))489PrintOut("%s", prompt);490}491492return FALSE;493}494495496void com_init(char *word, int flag_pw, int flag_ce, VARIABLE *(*sub)(),497int minp, int maxp, char *help_text )498/*======================================================================499? Adds commands to global command list.500|501| Parameters:502| char *word503| - the keyword user gives for this command to be executed.504| int flag_pw505| - flag telling if the command can be executed element506| by element using function *(*sub)().507| int flag_ce508| - flag telling if the command can be executed when509| preprocessing if constant arguments510| double *(*sub)()511| - function to be executed when this command is given512| int minp, maxp513| - maximum and minimum number of parameters to command514|515| The global list of available commands is updated (or created if516| nonexistent).517|518& lst_add()519~ *_com_init()520^=====================================================================*/521{522COMMAND *ptr; /* can't get rid of this */523524525/*526Fill the structure...527*/528ptr = (COMMAND *)ALLOCMEM(COMSIZE);529NAME(ptr) = STRCOPY(word);530if (flag_pw)531ptr->flags |= CMDFLAG_PW;532if (flag_ce)533ptr->flags |= CMDFLAG_CE;534ptr->minp = minp;535ptr->maxp = maxp;536ptr->sub = sub;537ptr->help = help_text;538539/*540...and update the list.541*/542lst_add(COMMANDS, (LIST *)ptr);543544return;545}546547void com_free(void)548/*======================================================================549? Deletes the list of commands and frees associated memory.550|551& lst_purge()552^=====================================================================*/553{554/*555Give memory back to system556*/557lst_purge(COMMANDS);558559return;560}561562COMMAND *com_check(char *str)563/*======================================================================564? Look for command from COMMANDS list by name.565|566= COMMAND *NULL if does not exist, pointer to command otherwise567& lst_find()568^=====================================================================*/569{570return (COMMAND *)lst_find(COMMANDS, str);571}572573VARIABLE *com_help( VARIABLE *ptr )574/*======================================================================575? Print list of commands and user defined functions from global lists.576|577! The command to get here is "help" but it really is not very helpful.578|579& lst_print()580^=====================================================================*/581{582COMMAND *cmd;583FUNCTION *fnc;584char *name;585586if ( !ptr )587{588589lst_print(COMMANDS);590lst_print(FUNCTIONS);591592} else {593594name = var_to_string( ptr );595596if ( (cmd = com_check( name ) ) != (COMMAND *)NULL )597{598599if ( cmd->help )600PrintOut( "\n%s\n", cmd->help );601else602PrintOut( "\nSorry: no help available on [%s].\n", name );603604} else if ( (fnc = fnc_check( name ) ) != (FUNCTION *)NULL )605{606607if ( fnc->help )608PrintOut( "\n%s", fnc->help );609else610PrintOut( "\nSorry: no help available on [%s].\n", name );611612} else {613614error( "help: symbol not found: [%s]\n", name );615616}617618FREEMEM( name );619}620621return (VARIABLE *)NULL;622}623624VARIABLE *com_pointw(double (*sub)(), VARIABLE *ptr)625/*======================================================================626? This routine does a function call (*sub)(), for each element in627| matrix given by ptr.628|629= a temporary VARIABLE for which M(res, i, j) = (*sub)(M(ptr, i, j)630& var_temp_new(), *(sub)()631^=====================================================================*/632{633VARIABLE *res,*ptr2; /* pointer to result structure */634635double *a, *a2, *a3, *b; /* pointer to matrices */636int n, m, sz; /* matrix dimensions */637638int i; /* loop index */639640/*641Get space for result and ...642*/643n = NROW(ptr); m = NCOL(ptr);644res = var_temp_new(TYPE(ptr) ,n , m);645646sz = n*m;647a = MATR(ptr); b = MATR(res);648649/*650...to action.651*/652ptr2 = NEXT(ptr);653if(ptr2)654{655if(n!=NROW(ptr2)||m!=NCOL(ptr2))656{657error("Pointwise function arguments must all be of same size.");658}659a2 = MATR(ptr2);660661ptr2 = NEXT(ptr2);662if(ptr2)663{664if(n!=NROW(ptr2)||m!=NCOL(ptr2))665{666error("Pointwise function arguments must all be of same size,");667}668if(NEXT(ptr2))669{670error("Currently at most three arguments for pointwise functions allowed, sorry.");671}672a3 = MATR(ptr2);673for(i = 0; i < sz; i++)674*b++ = ((double (*)(double, double, double)) sub)(*a++, *a2++, *a3++);675}676else677{678for(i = 0; i < sz; i++)679*b++ = ((double (*)(double, double)) sub)(*a++, *a2++);680}681}682else683{684for(i = 0; i < sz; i++)685*b++ = ((double (*)(double)) sub)(*a++);686}687688return res;689}690691VARIABLE *com_el(VARIABLE *ptr)692/*======================================================================693? Extracts specified elements from a matrix. Indexes are given by two694| column vectors. The values of the elements of these vectors give695| the required indexes. If there is only one index vector given696| it is assumed to be column index and row index is set to scalar 0.697|698| If matrix x is, for example,699|700| 1 2701| 3 4702|703| you get the first row by704|705| x[0, 0 1]706|707| or by708|709| x(0 1)710|711= A new temporary VARIABLE, whose size equals to712| number of row indexes times number of column indexes.713|714& var_temp_new(), var_delete_temp()715^=====================================================================*/716{717VARIABLE *res, /* result ... */718*par = NEXT(ptr); /* pointer to list of VARIABLES */719/* containing indexes */720721static double defind = 0.0;722#pragma omp threadprivate (defind)723double *ind1 = &defind, *ind2;724725int i, j, k, /* loop indexes */726rows, cols, /* no. of rows and columns in the matrix */727/* to be indexed. */728size1 = 1, size2,729ind;730731rows = NROW(ptr); cols = NCOL(ptr);732733/*734* check if scalar ....735*/736if (rows == 1 && cols == 1)737{738if (*MATR(par) != 0) error("Index out of bounds.\n");739if (NEXT(par) != NULL)740if (*MATR(NEXT(par)) != 0) error("Index out of bounds.\n");741res = var_temp_new(TYPE(ptr),1,1);742*MATR(res) = *MATR(ptr);743return res;744}745746/*747The matrix will be indexed by two column vectors.748If there is just one assume it's column index and749make rowindex 0.750*/751if (NEXT(par) == NULL)752{753if (NROW(par) == rows && NCOL(par) == cols)754{755int logical = TRUE,756onecount=0;757758double *dtmp;759760dtmp = MATR(par);761for(i = 0; i < NROW(par)*NCOL(par); i++)762if (dtmp[i] == 0)763{764}765else if (dtmp[i] == 1)766{767onecount++;768}769else770{771logical = FALSE;772break;773}774775if (logical)776{777if (onecount == 0) return NULL;778779res = var_temp_new(TYPE(ptr),1,onecount);780for(i=0,k=0; i < rows; i++)781for(j=0; j < cols; j++)782if (M(par,i,j) == 1)783{784memcpy(&M(res,0,k++),&M(ptr,i,j),sizeof(double));785}786return res;787}788}789790ind2 = MATR(par); size2 = NCOL(par);791cols *= rows; rows = 1;792}793else794{795ind1 = MATR(par); size1 = NCOL(par);796size2 = NCOL(NEXT(par));797ind2 = MATR(NEXT(par));798}799800/*801Space for result802*/803res = var_temp_new(TYPE(ptr), size1, size2);804805/*806Extract the values (try making sense out of that807if you feel like it).808*/809for(i = 0; i < size1; i++)810{811ind = (int)ind1[i];812for(j = 0; j < size2; j++)813if (ind < rows && (int)ind2[j] < cols)814memcpy(&M(res,i,j),&M(ptr,ind,(int)ind2[j]),sizeof(double));815else816error("Index out of bounds.\n");817}818819return res;820}821822VARIABLE *com_source(VARIABLE *ptr)823/*======================================================================824? Redirect input stream to a file, whose name is given.825|826@ FILE *math_in827& ALLOCMEM, FREEMEM, fopen(), fclose(), error()828^=====================================================================*/829{830char *name; /* Hold converted string (file name) */831832FILE *save_in = math_in; /* Save previous input stream until */833/* we are done with the new file. */834835/*836convert the file name from ptr.837*/838name = var_to_string(ptr);839840/*841Execute the file.842*/843if ((math_in = fopen(name,"r")) != NULL)844{845/* PrintOut("Executing commands from file, %s...\n", name); */846doread();847fclose(math_in);848}849else850{851PrintOut( "Source: Can't open file, %s.\n",name );852}853854math_in = save_in;855FREEMEM(name);856857return (VARIABLE *)NULL;858}859860861VARIABLE *com_apply(VARIABLE *ptr)862/*======================================================================863? Executes given string.864|865& ALLOCMEM, FREEMEM, doit()866^=====================================================================*/867{868VARIABLE *res; /* result pointer */869870char *p, *q; /* holds the string to be executed, after */871/* conversion from structure VARIABLE * */872873int i, j; /* just loop indexes */874875876/*877Allocate space for the string...878*/879p = q = ALLOCMEM(NROW(ptr) * NCOL(ptr) + 1);880881/*882... convert it ...883*/884for(i = 0; i < NROW(ptr); i++)885for(j = 0; j < NCOL(ptr); j++)886*p++ = (char)M(ptr,i,j);887888*p = '\0';889890/*891... and try executing it.892*/893res = doit( q );894895FREEMEM(q);896897return res;898}899900void mem_free(void *mem)901/*======================================================================902? Free memory given by argument, and unlink it from allocation list.903| Currently FREEMEM(ptr) is defined to be mem_free(ptr).904|905& free()906~ mem_alloc(), mem_free_all()907^=====================================================================*/908{909ALLOC_LIST *lst;910911#ifdef DEBUG912tot--; fprintf(fplog,"free addr: %d total: %d\n", ALLOC_LST(mem), tot);913fflush( fplog );914#endif915/*916if the list is empty return917*/918if ( (lst = (ALLOC_LIST *)ALLOC_HEAD) == (ALLOC_LIST *)NULL )919{920#if 1921/* ????? */922free( ALLOC_LST(mem) );923#else924fprintf( stderr, "SHOULD THIS HAPPEN ????\n" );925#endif926return;927}928929/*930* it's not the header, look if it's in list at all931*/932if (ALLOC_PTR(lst) != mem)933{934935for(; NEXT(lst); lst = NEXT(lst))936{937if (ALLOC_PTR(NEXT(lst)) == mem) break;938}939940/*941* item was not found from the list. free ptr and return.942*/943if (NEXT(lst) == (ALLOC_LIST *)NULL)944{945free(ALLOC_LST(mem));946return;947}948949/*950* unlink951*/952NEXT(lst) = NEXT(NEXT(lst));953}954955/*956* item was the header, unlink it957*/958else959ALLOC_HEAD = NEXT(ALLOC_HEAD);960961/*962* and at last return memory back to system963*/964free(ALLOC_LST(mem));965966return;967}968969void mem_free_all(void)970/*======================================================================971? Free all memory allocated since last entry of doread.972| (actually free all memory from list ALLOCATIONS).973|974~ mem_alloc(), mem_free(), doread(), error()975^=====================================================================*/976{977ALLOC_LIST *lst, *lstn;978979for(lst = (ALLOC_LIST *)ALLOC_HEAD; lst;)980{981#ifdef DEBUG982tot--; fprintf(fplog,"freeall addr: %d total: %d\n", lst, tot);983fflush( fplog );984#endif985lstn = NEXT(lst);986free( (char *)lst );987lst = lstn;988}989990ALLOC_HEAD = (LIST *)NULL; /* security */991992return;993}994995void *mem_alloc(size_t size)996/*======================================================================997? Allocate memory and link it to memory allocation list.998|999~ calloc(), free(), error()1000^=====================================================================*/1001{1002ALLOC_LIST *lst;10031004/*1005* try allocating memory1006*/1007if ((lst = (ALLOC_LIST *)calloc(size+sizeof(ALLOC_LIST), 1)) != NULL)1008{1009NEXT(lst) = (ALLOC_LIST *)ALLOC_HEAD; ALLOC_HEAD = (LIST *)lst;1010}1011else1012error("Can't alloc mem.\n");10131014#ifdef DEBUG1015tot++; fprintf(fplog,"alloc addr: %d size: %d total: %d\n",1016lst, size, tot);1017fflush( fplog );1018#endif1019return ALLOC_PTR(lst);1020}102110221023