Implementation of jacobian and invjacobian for the integrals

Author: softmattertheory

src/geometry/functional.c

@@ -4334,21 +4334,113 @@ static value integral_cgfn(vm *v, int nargs, value *args) {
     return out;
 }
 
+/* -------------------
+ * Jacobian
+ * ------------------- */
+
+/*
+ * A reference triangle is mapped to a target triangle through a
+ * linear transformation (the pushforward); an inverse transformation (the pullback)
+ * exists if the triangle is not degenerate. This function computes the forward
+ * and inverse jacobians.
+ */
+
+int jacobianhandle; // TL storage handle for Jacobian
+int invjacobianhandle; // TL storage handle for inverse Jacobian
+
+void _fetchvertices(objectintegralelementref *elref, objectmesh *mesh, int nv, elementid *vid, double **x) {
+    // Fetch reference vertices
+    for (int j=0; j<nv; j++) matrix_getcolumn(elref->iref->mref->vert, vid[j], &x[j]);
+}
+
+void _edgevectors(grade g, int dim, double **x, double *out) {
+    for (int i=0; i<g; i++) functional_vecsub(dim, x[i+1], x[0], out + i*dim);
+}
+
+/** Evaluates the jacobian and inverse jacobian; returns either of these as requested */
+void integral_evaluatejacobian(vm *v, value *jac, value *invjac) {
+    objectintegralelementref *elref = integral_getelementref(v);
+    
+    if (!elref) {
+        morpho_runtimeerror(v, INTEGRAL_SPCLFN, JACOBIAN_FUNCTION); return;
+    }
+    
+    int dim = elref->mesh->dim;     // Dimension of the mesh
+    
+    // Allocate matrices
+    objectmatrix *J=object_newmatrix(dim, dim, true);
+    objectmatrix *Jinv=object_newmatrix(dim, dim, true);
+    
+    if (J) vm_settlvar(v, jacobianhandle, MORPHO_OBJECT(J));
+    if (Jinv) vm_settlvar(v, invjacobianhandle, MORPHO_OBJECT(Jinv));
+    
+    if (!J || !Jinv) { morpho_runtimeerror(v, ERROR_ALLOCATIONFAILED); return; }
+    
+    // Now compute them
+    grade g = elref->g;             // Grade of the element
+    int nv = elref->nv;             //
+    
+    double **X = elref->vertexposn; // Vertex positions of the target element
+    double *x[nv];                  // Vertex positions of the reference element
+    
+    objectmesh *mref = elref->iref->mref; // Reference mesh
+    if (mref) _fetchvertices(elref, mref, nv, elref->vid, x);
+    
+    // Construct matrix of edge vectors for target and reference elements
+    double starget[dim*dim], sref[dim*dim], sinv[dim*dim];
+    objectmatrix St = MORPHO_STATICMATRIX(starget, dim, dim),
+                 Sr = MORPHO_STATICMATRIX(sref, dim, dim),
+                 Sinv = MORPHO_STATICMATRIX(sinv, dim, dim);
+    
+    _edgevectors(g, dim, X, starget);
+    if (mref) {
+        _edgevectors(g, dim, x, sref);
+        matrix_inverse(&Sr, &Sinv);
+    } else {
+        matrix_identity(&Sinv); // If no reference, the reference is the unit triangle
+    }
+    
+    matrix_mul(&St, &Sinv, J); // J = S . s^-1
+
+    matrix_inverse(J, Jinv); // Compute J^-1
+    
+    if (jac) *jac = MORPHO_OBJECT(J);
+    if (invjac) *invjac = MORPHO_OBJECT(Jinv);
+}
+
+static value integral_jacobian(vm *v, int nargs, value *args) {
+    value out=MORPHO_NIL;
+
+    vm_gettlvar(v, jacobianhandle, &out);
+    if (MORPHO_ISNIL(out)) integral_evaluatejacobian(v, &out, NULL);
+    
+    return out;
+}
+
+static value integral_invjacobian(vm *v, int nargs, value *args) {
+    value out=MORPHO_NIL;
+
+    vm_gettlvar(v, invjacobianhandle, &out);
+    if (MORPHO_ISNIL(out)) integral_evaluatejacobian(v, NULL, &out);
+    
+    return out;
+}
+
 /* ----------------------
  * General initialization
  * ---------------------- */
 
 /** Clears threadlocal storage */
 void integral_cleartlvars(vm *v) {
-    int handles[] = { elementhandle, normlhandle, tangenthandle, cauchygreenhandle, -1 };
+    int handles[] = { elementhandle, normlhandle, tangenthandle, cauchygreenhandle, jacobianhandle, invjacobianhandle, -1 };
 
     for (int i=0; handles[i]>=0; i++) {
         vm_settlvar(v, handles[i], MORPHO_NIL);
     }
 }
 
 void integral_freetlvars(vm *v) {
-    int handles[] = { normlhandle, tangenthandle, cauchygreenhandle, -1 };
+    int handles[] = { normlhandle, tangenthandle, cauchygreenhandle,jacobianhandle, invjacobianhandle, -1 };
 
     for (int i=0; handles[i]>=0; i++) {
         value val;
@@ -4952,6 +5044,8 @@ void functional_initialize(void) {
     builtin_addfunction(NORMAL_FUNCTION, integral_normal, BUILTIN_FLAGSEMPTY);
     builtin_addfunction(GRAD_FUNCTION, integral_gradfn, BUILTIN_FLAGSEMPTY);
     builtin_addfunction(CGTENSOR_FUNCTION, integral_cgfn, BUILTIN_FLAGSEMPTY);
+    builtin_addfunction(JACOBIAN_FUNCTION, integral_jacobian, BUILTIN_FLAGSEMPTY);
+    builtin_addfunction(INVJACOBIAN_FUNCTION, integral_invjacobian, BUILTIN_FLAGSEMPTY);
 
     morpho_defineerror(VOLUMEENCLOSED_ZERO, ERROR_HALT, VOLUMEENCLOSED_ZERO_MSG);
     morpho_defineerror(FUNC_INTEGRAND_MESH, ERROR_HALT, FUNC_INTEGRAND_MESH_MSG);
@@ -4990,6 +5084,8 @@ void functional_initialize(void) {
     tangenthandle=vm_addtlvar();
     normlhandle=vm_addtlvar();
     cauchygreenhandle=vm_addtlvar();
+    jacobianhandle=vm_addtlvar();
+    invjacobianhandle=vm_addtlvar();
 
     morpho_addfinalizefn(functional_finalize);
 }

src/geometry/functional.h

@@ -60,6 +60,8 @@
 #define NORMAL_FUNCTION                "normal"
 #define GRAD_FUNCTION                  "grad"
 #define CGTENSOR_FUNCTION              "cgtensor"
+#define JACOBIAN_FUNCTION              "jacobian"
+#define INVJACOBIAN_FUNCTION           "invjacobian"
 
 /* Functional names */
 #define LENGTH_CLASSNAME               "Length"

test/functionals/areaintegral/jacobian.morpho

@@ -0,0 +1,33 @@
+// Test the jacobian() function 
+
+import meshtools 
+
+// Create the reference triangle
+var mb = MeshBuilder() 
+mb.addvertex([0,0])
+mb.addvertex([1,0])
+mb.addvertex([0,1])
+mb.addedge([0,1])
+mb.addedge([1,2])
+mb.addedge([2,0])
+mb.addface([0,1,2])
+
+var mref = mb.build()
+
+// Make the target triangle by cloning the reference..
+var m = mref.clone() 
+m.setvertexmatrix(2*m.vertexmatrix()) // ... and scaling by 2
+
+fn integrand(x) {
+    var F = jacobian() 
+    //var Jinv = invjacobian()
+    var Ft = F.transpose() 
+
+    return (F*Ft).trace()
+}
+
+var a=AreaIntegral(integrand, reference=mref, weightByReference=true).total(m)
+
+print a
+// expect: 4
+// because Tr(F*Ft)=8; and the area of the reference tri. is 1/2