Correct electro/magnetostatic force for axisymmetry

examples/cmut-collapse-mode-axisymmetry-2d/main.cpp

@@ -109,7 +109,7 @@ void sparselizard(void)
     // the gradient of the previously computed electric potential field and the electric permittivity.
     //
     // The electrostatic forces are computed on the mesh deformed by field umesh.
-    elasticity += integral(wholedomain, umesh, predefinedelectrostaticforce(grad(tf(u,solid)), grad(v), epsilon));
+    elasticity += integral(wholedomain, umesh, predefinedelectrostaticforce(tf(u,solid), grad(v), epsilon));
 
 
     // Solve the Laplace equation in the cavity to smoothly deform the mesh.
@@ -209,7 +209,7 @@ void sparselizard(void)
     helasticity += integral(solid, predefinedelasticity(dof(uh), tf(uh), u, E, nu, 0.0));
     // Here the electrostatic force must be computed using an FFT (on 5 timesteps)
     // because the force calculation involves nonlinear operations on multiharmonic fields:
-    helasticity += integral(wholedomain, 5, umesh, predefinedelectrostaticforce(grad(tf(uh,solid)), grad(vh), epsilon));
+    helasticity += integral(wholedomain, 5, umesh, predefinedelectrostaticforce(tf(uh,solid), grad(vh), epsilon));
     // The inertia term is added for the harmonic analysis:
     helasticity += integral(solid, -rho*dtdt(dof(uh))*tf(uh));
 
@@ -231,7 +231,7 @@ void sparselizard(void)
     std::cout << "Peak AC deflection: " << uacmax*1e9 << " nm" << std::endl;
 
     // Code validation line. Can be removed.
-    std::cout << (uacmax < 0.977925e-9 && uacmax > 0.977923e-9);
+    std::cout << (uacmax < 0.977678e-9 && uacmax > 0.977675e-9);
 }
 
 

examples/magnetostatics-scalar-potential-2d/main.cpp

@@ -116,7 +116,7 @@ void sparselizard(void)
 	// The force term integral must be performed on all elements in the steel region AND in the ELEMENT LAYER AROUND.
 	// This is obtained by considering the elements on the whole domain 'wholedomain'. Since the force must only be
 	// calculated on the steel disk, the test functions are only selected on the steel (this is done with tf(..., steel)). 
-	forceprojection += integral(wholedomain, - predefinedmagnetostaticforce(grad(tf(magforce, steel)), -grad(phi), mu));
+	forceprojection += integral(wholedomain, - predefinedmagnetostaticforce(tf(magforce, steel), -grad(phi), mu));
 
 	forceprojection.generate();
 	vec solf = solve(forceprojection.A(), forceprojection.b());

examples/nonlinear-electro-elasto-acoustic-micromembrane-multiharmonic-2d/main.cpp

@@ -139,7 +139,8 @@ void sparselizard(void)
     //
     // The electrostatic forces must be computed on the mesh deformed by field umesh.
     // The calculations are brought back to the undeformed configuration.
-    elastoacoustic += integral(electricdomain, 20, predefinedelectrostaticforce(transpose(invjac*transpose(grad(tf(u,solid)))), invjac*grad(v), epsilon) * detjac );
+    expression gradtfu = invjac*expression({{grad(compx(tf(u,solid)))}, {grad(compy(tf(u,solid)))}});
+    elastoacoustic += integral(electricdomain, 20, predefinedelectrostaticforce({gradtfu.getcolumn(0),gradtfu.getcolumn(1)}, invjac*grad(v), epsilon) * detjac );
 
     // The wave equation is solved in the fluid:
     elastoacoustic += integral(fluid, -grad(dof(p))*grad(tf(p)) -1/pow(c,2)*dtdt(dof(p))*tf(p));

examples/nonlinear-electroelastic-micromembrane-2d/main.cpp

@@ -89,7 +89,7 @@ void sparselizard(void)
     // the gradient of the previously computed electric potential field and the electric permittivity.
     //
     // The electrostatic forces are computed on the mesh deformed by field umesh.
-    elasticity += integral(electricdomain, umesh, predefinedelectrostaticforce(grad(tf(u,solid)), grad(v), epsilon));
+    elasticity += integral(electricdomain, umesh, predefinedelectrostaticforce(tf(u,solid), grad(v), epsilon));
 
 
     // Solve the Laplace equation in the vacuum gap to smoothly deform the mesh.

src/expression/operation/mathop.cpp

@@ -790,26 +790,52 @@ expression mathop::predefinedelasticity(expression dofu, expression tfu, field u
     }
 }
 
-expression mathop::predefinedelectrostaticforce(expression gradtfu, expression E, expression epsilon)
+expression mathop::predefinedelectrostaticforce(expression tfu, expression E, expression epsilon)
+{
+    if (tfu.countcolumns() != 1)
+    {
+        std::cout << "Error in 'mathop' namespace: the force formula expected a column vector expression as first argument" << std::endl;
+        abort();
+    }
+
+	std::vector<expression> spacederivatives(tfu.countrows());
+	for (int i = 0; i < tfu.countrows(); i++)
+		spacederivatives[i] = grad(tfu.at(i,0));
+
+    return predefinedelectrostaticforce(spacederivatives, E, epsilon);
+}
+
+expression mathop::predefinedelectrostaticforce(std::vector<expression> dxyztfu, expression E, expression epsilon)
 {
     E.reuseit();
+    epsilon.reuseit();
+
+    std::vector<std::vector<expression>> exprs(dxyztfu.size());
+    for (int i = 0; i < dxyztfu.size(); i++)
+    	exprs[i] = {dxyztfu[i]};
+
+    // Scalar gradient here:
+    expression gradtfu(exprs);
 
     if (gradtfu.countcolumns() == 1)
     {
-        std::cout << "Error in 'mathop' namespace: force formula is undefined for 1D displacements" << std::endl;
+        std::cout << "Error in 'mathop' namespace: the force formula is undefined for 1D displacements" << std::endl;
         abort();
     }
 
-    gradtfu = gradtfu.transpose();
-
     if (gradtfu.countcolumns() == 2)
         return -( epsilon*0.5 * (pow(compx(E),2) * entry(0,0,gradtfu) - pow(compy(E),2) * entry(0,0,gradtfu) + 2 * compx(E) * compy(E) * entry(1,0,gradtfu))      +epsilon*0.5 * (-pow(compx(E),2) * entry(1,1,gradtfu) + pow(compy(E),2) * entry(1,1,gradtfu) + 2 * compy(E) * compx(E) * entry(0,1,gradtfu)) );
     if (gradtfu.countcolumns() == 3)
         return -( epsilon*0.5 * (pow(compx(E),2) * entry(0,0,gradtfu) - pow(compy(E),2) * entry(0,0,gradtfu) - pow(compz(E),2) * entry(0,0,gradtfu) + 2 * compx(E) * compy(E) * entry(1,0,gradtfu) + 2 * compx(E) * compz(E) * entry(2,0,gradtfu))      +epsilon*0.5 * (-pow(compx(E),2) * entry(1,1,gradtfu) + pow(compy(E),2) * entry(1,1,gradtfu) - pow(compz(E),2) * entry(1,1,gradtfu) + 2 * compy(E) * compx(E) * entry(0,1,gradtfu) + 2 * compy(E) * compz(E) * entry(2,1,gradtfu))      +epsilon*0.5 * (-pow(compx(E),2) * entry(2,2,gradtfu) - pow(compy(E),2) * entry(2,2,gradtfu) + pow(compz(E),2) * entry(2,2,gradtfu) + 2 * compz(E) * compx(E) * entry(0,2,gradtfu) + 2 * compz(E) * compy(E) * entry(1,2,gradtfu)) );
 }
 
-expression mathop::predefinedmagnetostaticforce(expression gradtfu, expression H, expression mu)
+expression mathop::predefinedmagnetostaticforce(expression tfu, expression H, expression mu)
+{
+	return predefinedelectrostaticforce(tfu, H, mu);
+}
+
+expression mathop::predefinedmagnetostaticforce(std::vector<expression> dxyztfu, expression H, expression mu)
 {
-	return predefinedelectrostaticforce(gradtfu, H, mu);
+	return predefinedelectrostaticforce(dxyztfu, H, mu);
 }
 

src/expression/operation/mathop.h

@@ -151,8 +151,10 @@ namespace mathop
     // General anisotropic elasticity with geometrical nonlinearity and prestress (ignored if zero):
     expression predefinedelasticity(expression dofu, expression tfu, field u, expression hookesmatrix, expression prestress, std::string myoption = "");
 
-    expression predefinedelectrostaticforce(expression gradtfu, expression E, expression epsilon);
-    expression predefinedmagnetostaticforce(expression gradtfu, expression H, expression mu);
+    expression predefinedelectrostaticforce(expression tfu, expression E, expression epsilon);
+    expression predefinedelectrostaticforce(std::vector<expression> dxyztfu, expression E, expression epsilon);
+    expression predefinedmagnetostaticforce(expression tfu, expression H, expression mu);
+    expression predefinedmagnetostaticforce(std::vector<expression> dxyztfu, expression H, expression mu);
 };
 
 #endif