READ_DATA

Author: gtheler

ChangeLog.md

@@ -1,6 +1,7 @@
 # Current (under development)
 
  * ...
+ * keyword `READ_DATA` to read variables and vectors from files
  * `PROBLEM` can define `MESH` in a single line without an explicit `READ_MESH` instruction
  * Neo-hookean material model
  * MMS verification for mechanical with both sdef & ldef

Makefile.am

@@ -47,6 +47,7 @@ TESTS_general = \
         tests/petsc_options.sh \
         tests/point-kinetics.sh \
         tests/print.sh \
+        tests/read.sh \
         tests/solve.sh \
         tests/trig.sh \
         tests/transient-mesh.sh \

src/Makefile-base.am

@@ -50,6 +50,7 @@ feenox_SOURCES = \
 ./parser/parser.h \
 ./parser/auxiliary.c \
 ./io/file.c \
+./io/read.c \
 ./io/print.c \
 ./io/printf.c \
 ./math/alias.c \

src/feenox.h

@@ -329,6 +329,8 @@ typedef struct printf_t printf_t;
 typedef struct print_function_t print_function_t;
 typedef struct multidim_range_t multidim_range_t;
 typedef struct print_vector_t print_vector_t;
+typedef struct read_data_t read_data_t;
+
 
 typedef struct sort_vector_t sort_vector_t;
 typedef struct phase_object_t phase_object_t;
@@ -830,6 +832,15 @@ struct printf_t {
   printf_t *next;
 };
 
+struct read_data_t {
+  file_t *file;
+  // TODO: several things mixed together, this is where I miss C++'s STL
+  var_t *variable;
+  vector_t *vector;
+  
+  read_data_t *next;
+};
+
 struct multidim_range_t {
   int dimensions;
 
@@ -2088,6 +2099,10 @@ FILE *feenox_fopen(const char *filepath, const char *mode);
 extern int feenox_instruction_file_open(void *arg);
 extern int feenox_instruction_file_close(void *arg);
 
+// read.c
+extern int feenox_instruction_read_data(void *arg);
+
+
 // abort.c
 extern int feenox_instruction_abort(void *arg);
 

src/io/read.c

@@ -0,0 +1,42 @@
+/*------------ -------------- -------- --- ----- ---   --       -            -
+ *  feenox read data from files
+ *
+ *  Copyright (C) 2025 Jeremy Theler
+ *
+ *  This file is part of feenox.
+ *
+ *  feenox is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  FeenoX is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with feenox.  If not, see <http://www.gnu.org/licenses/>.
+ *------------------- ------------  ----    --------  --     -       -         -
+ */
+#include "feenox.h"
+
+int feenox_instruction_read_data(void *arg) {
+ 
+  read_data_t *read_data = (read_data_t *)arg;
+  if (read_data->file->pointer == NULL) {
+    feenox_call(feenox_instruction_file_open(read_data->file));
+  }  
+  
+  if (read_data->variable != NULL) {
+    fscanf(read_data->file->pointer, "%lf", read_data->variable->value);
+  } else {
+    double xi = 0;
+    for (int i = 0; i < read_data->vector->size; i++) {
+      fscanf(read_data->file->pointer, "%lf", &xi);
+      gsl_vector_set(read_data->vector->value, i, xi);
+    }
+  }
+  
+  return FEENOX_OK;
+}

src/math/function.c

@@ -513,6 +513,19 @@ double feenox_function_eval(function_t *this, const double *const_x) {
 
   // check if we need to initialize
   if (this->initialized == 0) {
+    // there is a pathological case where a BC or a property depends on T(x,y,z)
+    // but T(x,y,z) itself is not ready before SOLVE_PROBLEM so we have to look
+    // for the initial condition T_0(x,y,z)...
+    if (this->algebraic_expression.items == NULL && this->mesh == NULL && this->vector_argument == NULL) {
+      function_t *initial = feenox_get_function_ptr("T_0");
+      if (initial != NULL) {
+        return feenox_function_eval(initial, x);
+      } else {
+        feenox_push_error_message("function '%s' is not ready to be used and no initial condition '%s_0' found", this->name, this->name);
+        return FEENOX_ERROR;
+      }
+    }
+      
     if (feenox_function_init(this) != FEENOX_OK) {
       feenox_runtime_error();
     }

src/parser/parser.c

@@ -191,7 +191,13 @@ int feenox_parse_line(void) {
     } else if (strcasecmp(token, "FUNCTION") == 0) {
       feenox_call(feenox_parse_function());
       return FEENOX_OK;
-    
+
+///kw+READ_DATA+desc Read data from a file (or standard output, or pipe, or queue, etc.).
+///kw+PRINT+usage READ_DATA
+    } else if (strcasecmp(token, "READ_DATA") == 0) {
+      feenox_call(feenox_parse_read_data());
+      return FEENOX_OK;
+      
 ///kw+PRINT+desc Write plain-text and/or formatted data to the standard output or into an output file.
 ///kw+PRINT+usage PRINT
     } else if (strcasecmp(token, "PRINT") == 0) {
@@ -1568,6 +1574,34 @@ int feenox_parse_open_close(const char *what) {
   return FEENOX_OK;  
 }
 
+int feenox_parse_read_data(void) {
+
+  read_data_t *read_data = NULL;
+  feenox_check_alloc(read_data = calloc(1, sizeof(read_data_t)));
+  
+  char *token = NULL;
+  while ((token = feenox_get_next_token(NULL)) != NULL) {
+///kw+READ_DATA+usage [ FILE { <file_path> | <file_id> } ]
+    if (strcasecmp(token, "FILE") == 0 || strcasecmp(token, "FILE_PATH") == 0) {
+      feenox_call(feenox_parser_file(&read_data->file));
+      if (read_data->file->mode == NULL) {
+        feenox_check_alloc(read_data->file->mode = strdup("r"));
+      }
+    } else {
+      if ((read_data->variable = feenox_get_variable_ptr(token)) == NULL) {
+        if ((read_data->vector = feenox_get_vector_ptr(token)) == NULL) {
+          feenox_push_error_message("unknown symbol '%s'", token);
+          return FEENOX_ERROR;
+        }
+      }
+      
+    }
+  }
+  
+  feenox_call(feenox_add_instruction(feenox_instruction_read_data, read_data));
+  
+  return FEENOX_OK;
+}
 int feenox_parse_print(void) {
 
   // I don't expect anybody to want to use this PRINT instruction through the API

src/parser/parser.h

@@ -100,6 +100,8 @@ extern int feenox_parse_sort_vector(void);
 extern int feenox_parse_file(char *mode);
 extern int feenox_parse_open_close(const char *what);
 
+extern int feenox_parse_read_data(void);
+
 extern int feenox_parse_print(void);
 extern int feenox_parse_printf(int all_ranks);
 extern int feenox_parse_print_function(void);

src/pdes/thermal/bc.c

@@ -64,6 +64,7 @@ int feenox_problem_bc_parse_thermal(bc_data_t *bc_data, const char *lhs, char *r
 
   // for non-linear problems it is important to have a good initial guess
   // if the user did not give us one in T_guess we average all the temperatures from the BCs
+  // but we have to be sure that the BC expression does not involve a function coming from somewhere else
   if (bc_data->type_phys == BC_TYPE_THERMAL_TEMPERATURE || bc_data->type_phys == BC_TYPE_THERMAL_CONVECTION_TEMPERATURE ) {
     thermal.guessed_initial_guess += feenox_expression_eval(&bc_data->expr);
     thermal.n_bc_temperatures++;

tests/onetwothree.txt

@@ -0,0 +1 @@
+123 1 2 3