Merge pull request #234 from Morpho-lang/mapfieldgradient

functional_mapfieldgradient

examples/examples.py

@@ -86,7 +86,7 @@ def run(file,testLog,CI):
 
             if (iserror(line)):
                 if not CI:
-                    print(stylize("Failed",colored.fg("red")))
+                    print("Failed") #stylize("Failed",colored.fg("red"))) // Temporarily disable this 6/19/23 due to colored module API change
                 else:
                     print("::error file = {",file,"}::{",file," Failed}")
 
@@ -101,7 +101,7 @@ def run(file,testLog,CI):
         if (ret ==1):
             if not CI:
                 print(file+":", end=" ")
-                print(stylize("Passed",colored.fg("green")))
+                print("Passed") #stylize("Passed",colored.fg("green")))
         # Delete the temporary file
         os.system('rm ' + tmp)
     return ret

morpho5/Makefile

@@ -50,7 +50,7 @@ headers:
 	find . -type f -name '*.h' -exec cp {} $(HEADERSDIR) \;
 
 .PHONY: extensions
-headers: 
+extensions: 
 	mkdir -p $(EXTENSIONSDIR)
 
 .PHONY: install

morpho5/Makefile.linux

@@ -50,7 +50,7 @@ headers:
 	find . -type f -name '*.h' -exec cp {} $(HEADERSDIR) \;
 
 .PHONY: extensions
-headers: 
+extensions: 
 	mkdir -p $(EXTENSIONSDIR)
 
 .PHONY: install

morpho5/Makefile.m1

@@ -50,7 +50,7 @@ headers:
 	find . -type f -name '*.h' -exec cp {} $(HEADERSDIR) \;
 
 .PHONY: extensions
-headers: 
+extensions: 
 	mkdir -p $(EXTENSIONSDIR)
 
 .PHONY: install

morpho5/geometry/functional.c

@@ -380,6 +380,72 @@ bool functional_mapgradientX(vm *v, functional_mapinfo *info, value *out) {
     return ret;
 }
 
+/** Calculate field gradient
+ * @param[in] v - virtual machine in use
+ * @param[in] info - map info structure
+ * @param[out] out - a field of integrand values
+ * @returns true on success, false otherwise. Error reporting through VM. */
+bool functional_mapfieldgradientX(vm *v, functional_mapinfo *info, value *out) {
+    objectmesh *mesh = info->mesh;
+    objectfield *field = info->field;
+    objectselection *sel = info->sel;
+    objectfield *grad=NULL;
+    grade g = info->g;
+    functional_fieldgradient *fgrad = info->fieldgrad;
+    void *ref = info->ref;
+    objectsparse *s=NULL;
+    bool ret=false;
+    int n=0;
+
+    /* How many elements? */
+    if (!functional_countelements(v, mesh, g, &n, &s)) return false;
+
+    /* Create the output field */
+    if (n>0) {
+        grad=object_newfield(mesh, field->prototype, field->dof);
+        if (!grad) { morpho_runtimeerror(v, ERROR_ALLOCATIONFAILED); return false; }
+        field_zero(grad);
+    }
+
+    if (grad) {
+        int vertexid; // Use this if looping over grade 0
+        int *vid=(g==0 ? &vertexid : NULL),
+            nv=(g==0 ? 1 : 0); // The vertex indices
+
+
+        if (sel) { // Loop over selection
+            if (sel->selected[g].count>0) for (unsigned int k=0; k<sel->selected[g].capacity; k++) {
+                if (!MORPHO_ISINTEGER(sel->selected[g].contents[k].key)) continue;
+
+                elementid i = MORPHO_GETINTEGERVALUE(sel->selected[g].contents[k].key);
+                if (s) sparseccs_getrowindices(&s->ccs, i, &nv, &vid);
+                else vertexid=i;
+
+                if (vid && nv>0) {
+                    if (!(*fgrad) (v, mesh, i, nv, vid, ref, grad)) goto functional_mapfieldgradient_cleanup;
+                }
+            }
+        } else { // Loop over elements
+            for (elementid i=0; i<n; i++) {
+                if (s) sparseccs_getrowindices(&s->ccs, i, &nv, &vid);
+                else vertexid=i;
+
+                if (vid && nv>0) {
+                    if (!(*fgrad) (v, mesh, i, nv, vid, ref, grad)) goto functional_mapfieldgradient_cleanup;
+                }
+            }
+        }
+
+        *out = MORPHO_OBJECT(grad);
+        ret=true;
+    }
+
+functional_mapfieldgradient_cleanup:
+    if (!ret) object_free((object *) grad);
+
+    return ret;
+}
+
 static bool functional_numericalremotegradient(vm *v, functional_mapinfo *info, objectsparse *conn, elementid remoteid, elementid i, int nv, int *vid, objectmatrix *frc);
 
 /* Calculates a numerical gradient */
@@ -1277,6 +1343,57 @@ bool functional_mapnumericalgradient(vm *v, functional_mapinfo *info, value *out
     return success;
 }
 
+/* ----------------------------
+ * Map field gradients
+ * ---------------------------- */
+
+/** Compute the field gradient */
+bool functional_mapfieldgradient(vm *v, functional_mapinfo *info, value *out) {
+    int success=false;
+    int ntask=morpho_threadnumber();
+    if (!ntask) return functional_mapfieldgradientX(v, info, out);
+    functional_task task[ntask];
+
+    varray_elementid imageids;
+    varray_elementidinit(&imageids);
+
+    objectfield *new[ntask];
+    for (int i=0; i<ntask; i++) new[i]=NULL;
+
+    if (!functional_preparetasks(v, info, ntask, task, &imageids)) return false;
+
+    /* Create output fields */
+    for (int i=0; i<ntask; i++) {
+        // Create one per thread
+        new[i]=object_newfield(info->mesh, info->field->prototype, info->field->dof);
+        if (!new[i]) { morpho_runtimeerror(v, ERROR_ALLOCATIONFAILED); goto functional_mapfieldgradient_cleanup; }
+        field_zero(new[i]);
+
+        task[i].mapfn=(functional_mapfn *) info->fieldgrad;
+        task[i].result=(void *) new[i];
+    }
+
+    functional_parallelmap(ntask, task);
+
+    /* Then add up all the fields using their underlying data stores */
+    for (int i=1; i<ntask; i++) matrix_add(&new[0]->data, &new[1]->data, &new[0]->data);
+
+    // TODO: Use symmetry actions
+    //if (info->sym==SYMMETRY_ADD) functional_symmetrysumforces(info->mesh, new[0]);
+
+    success=true;
+
+functional_mapfieldgradient_cleanup:
+    for (int i=1; i<ntask; i++) if (new[i]) object_free((object *) new[i]);
+
+    // ...and return the result
+    *out = MORPHO_OBJECT(new[0]);
+
+    functional_cleanuptasks(v, ntask, task);
+    varray_elementidclear(&imageids);
+    return success;
+}
+
 /* ----------------------------
  * Map numerical field gradients
  * ---------------------------- */
@@ -3081,7 +3198,7 @@ bool gradsq_evaluategradient(objectmesh *mesh, objectfield *field, int nv, int *
  @param[in] vid - vertex ids
  @param[out] out - should be field->psize * mesh->dim units of storage */
 bool gradsq_evaluategradient1d(objectmesh *mesh, objectfield *field, int nv, int *vid, double *out) {
-    UNREACHABLE("GradSq in 3D not implemented.");
+    UNREACHABLE("GradSq in 1D not implemented.");
     double *f[nv]; // Field value lists
     double *x[nv]; // Vertex coordinates
     unsigned int nentries=0;
@@ -3250,6 +3367,7 @@ value GradSq_fieldgradient(vm *v, int nargs, value *args) {
             info.cloneref = gradsq_cloneref;
             info.ref = &ref;
             functional_mapnumericalfieldgradient(v, &info, &out);
+            //functional_mapfieldgradient(v, &info, &out);
         } else morpho_runtimeerror(v, GRADSQ_ARGS);
     }
     if (!MORPHO_ISNIL(out)) morpho_bindobjects(v, 1, &out);

morpho5/geometry/functional.h

@@ -154,6 +154,9 @@ typedef bool (functional_integrand) (vm *v, objectmesh *mesh, elementid id, int
 /** Gradient function */
 typedef bool (functional_gradient) (vm *v, objectmesh *mesh, elementid id, int nv, int *vid, void *ref, objectmatrix *frc);
 
+/** Field gradient function */
+typedef bool (functional_fieldgradient) (vm *v, objectmesh *mesh, elementid id, int nv, int *vid, void *ref, objectfield *frc);
+
 struct s_functional_mapinfo; // Resolve circular typedef dependency
 
 /** Clone reference function */
@@ -173,6 +176,7 @@ typedef struct s_functional_mapinfo {
     elementid id; // Element id at which to evaluate the integrand
     functional_integrand *integrand; // Integrand function
     functional_gradient *grad; // Gradient
+    functional_fieldgradient *fieldgrad; // Field gradient
     functional_dependencies *dependencies; // Dependencies
     functional_cloneref *cloneref; // Clone a reference with a given field substituted
     functional_freeref *freeref; // Free a reference
@@ -190,6 +194,7 @@ bool functional_sumintegrand(vm *v, functional_mapinfo *info, value *out);
 bool functional_mapintegrand(vm *v, functional_mapinfo *info, value *out);
 bool functional_mapintegrandat(vm *v, functional_mapinfo *info, value *out);
 bool functional_mapgradient(vm *v, functional_mapinfo *info, value *out);
+bool functional_mapfieldgradient(vm *v, functional_mapinfo *info, value *out);
 bool functional_mapnumericalgradient(vm *v, functional_mapinfo *info, value *out);
 bool functional_mapnumericalfieldgradient(vm *v, functional_mapinfo *info, value *out);