diff --git a/CircuitCoupling/R_circuit.geo b/CircuitCoupling/R_circuit.geo
new file mode 100644
index 0000000000000000000000000000000000000000..8a5b6f29daa0dbfc9fcddfe18494f921cc38f210
--- /dev/null
+++ b/CircuitCoupling/R_circuit.geo
@@ -0,0 +1,39 @@
+// Just 4 cubes
+
+Include "R_circuit_common.pro";
+
+lc = 0.25; // Cjaracteristic length
+
+Point(1) = {0, 0, 0, lc};
+Extrude {1, 0, 0} {
+ Point{1};
+}
+Extrude {0, 1, 0} {
+ Line{1}; // Layers{10}; Recombine;
+}
+Extrude {0, 0, 1} {
+ Surface{5}; // Layers{10}; Recombine;
+}
+Translate {2, 0, 0} {
+ Duplicata { Volume{1}; }
+}
+Translate {0, 4, 0} {
+ Duplicata { Volume{1}; }
+}
+Translate {2, 4, 0} {
+ Duplicata { Volume{1}; }
+}
+
+Physical Surface(Cube1Top) = {27};
+Physical Surface(Cube1Bottom) = {5};
+Physical Surface(Cube2Top) = {34};
+Physical Surface(Cube2Bottom) = {29};
+Physical Surface(Cube3Top) = {61};
+Physical Surface(Cube3Bottom) = {56};
+Physical Surface(Cube4Top) = {88};
+Physical Surface(Cube4Bottom) = {83};
+
+Physical Volume(Cube1) = {1};
+Physical Volume(Cube2) = {28};
+Physical Volume(Cube3) = {55};
+Physical Volume(Cube4) = {82};
diff --git a/CircuitCoupling/R_circuit.jpg b/CircuitCoupling/R_circuit.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..14798d1d1a5249c1644f19b313a5cf35fcd967a7
Binary files /dev/null and b/CircuitCoupling/R_circuit.jpg differ
diff --git a/CircuitCoupling/R_circuit.pro b/CircuitCoupling/R_circuit.pro
new file mode 100644
index 0000000000000000000000000000000000000000..c93b64254e5510249146ae6b6c77bb473652bdfa
--- /dev/null
+++ b/CircuitCoupling/R_circuit.pro
@@ -0,0 +1,334 @@
+/* -------------------------------------------------------------------
+ Tutorial 8a : circuit coupling - resistive circuit
+
+ Features:
+ - Electrokinetic formulation coupled with resistive circuit
+ - Definition of circuit elements (sources, resistances)
+ - Implementation of a netlist
+
+ To compute the solution in a terminal:
+ getdp R_circuit -solve Jstatic -pos Cube_Top_Values_as_ASCII_File
+
+ To compute the solution interactively from the Gmsh GUI:
+ File > Open > R_circuit.pro
+ Run (button at the bottom of the left panel)
+ ------------------------------------------------------------------- */
+
+Include "R_circuit_common.pro";
+
+// FEM domain group data
+Group {
+ // Surfaces
+ Top1 = Region[Cube1Top];
+ Bottom1 = Region[Cube1Bottom];
+ Top2 = Region[Cube2Top];
+ Bottom2 = Region[Cube2Bottom];
+ Top3 = Region[Cube3Top];
+ Bottom3 = Region[Cube3Bottom];
+ Top4 = Region[Cube4Top];
+ Bottom4 = Region[Cube4Bottom];
+
+ // Volumes
+ Cube1 = Region[Cube1];
+ Cube2 = Region[Cube2];
+ Cube3 = Region[Cube3];
+ Cube4 = Region[Cube4];
+
+ // FEM Electrical regions
+ // * all volumes:
+ Vol_Ele = Region[{Cube1, Cube2, Cube3, Cube4}];
+ // * all surfaces connected to the lumped element circuit
+ Sur_Ele = Region[{Top1, Bottom1, Top2, Bottom2, Top3, Bottom3, Top4, Bottom4}];
+ // * total electrical computation domain:
+ VolWithSur_Ele = Region[{Vol_Ele, Sur_Ele}];
+}
+
+// Circuit group data (fictive regions for the circuit elements)
+Group {
+ // Sources
+ CurrentSource1 = Region[{3001}];
+ VoltageSource1 = Region[{3002}];
+
+ SourceI_Cir = Region[{CurrentSource1}]; // all current sources
+ SourceV_Cir = Region[{VoltageSource1}]; // all voltage sources
+ Sources_Cir = Region[{SourceI_Cir, SourceV_Cir}];
+
+ // Resistors
+ Res1 = Region[{4001}];
+ Res2 = Region[{4002}];
+ Res3 = Region[{4003}];
+ Res4 = Region[{4004}];
+
+ Resistance_Cir = Region[{Res1, Res2, Res3, Res4}]; // All resistors
+
+ // Complete circuit domain containing all circuit elements
+ Domain_Cir = Region[{Resistance_Cir, SourceI_Cir, SourceV_Cir}];
+}
+
+// FEM domain function data
+Function {
+ // Absolute voltage at the bottom of all cubes; set via constraint Voltage_3D,
+ // not via a voltage source
+ Vbottom = 4.0;
+
+ // Geometry of a cube
+ w = 1.0; // Width [m]
+ l = 1.0; // Length [m]
+ h = 1.0; // Height [m]
+
+ // Resistance
+ Rcube1 = 20.0; // Resistance of cube 1 [Ohm]
+ Rcube2 = 40.0; // Resistance of cube 2 [Ohm]
+ Rcube3 = 60.0; // Resistance of cube 3 [Ohm]
+ Rcube4 = 10.0; // Resistance of cube 4 [Ohm]
+
+ // Specific electrical conductivity [S*m/m^2]
+ kappa[Cube1] = h / (w * l * Rcube1);
+ kappa[Cube2] = h / (w * l * Rcube2);
+ kappa[Cube3] = h / (w * l * Rcube3);
+ kappa[Cube4] = h / (w * l * Rcube4);
+}
+
+// Circuit domain function data
+Function {
+ // External impedances functions
+ R[Res1] = 10.0;
+ R[Res2] = 20.0;
+ R[Res3] = 15.0;
+ R[Res4] = 10.0;
+}
+
+// FEM domain constraint data
+Constraint {
+ { Name Current_3D ; Type Assign ;
+ Case {
+ }
+ }
+ { Name Voltage_3D ; Type Assign ;
+ Case {
+ { Region Bottom1; Type Assign; Value Vbottom; }
+ { Region Bottom2; Type Assign; Value Vbottom; }
+ { Region Bottom3; Type Assign; Value Vbottom; }
+ { Region Bottom4; Type Assign; Value Vbottom; }
+ }
+ }
+}
+
+// Circuit domain constraint data
+Constraint {
+ { Name Current_Cir ;
+ Case {
+ { Region CurrentSource1; Value 1.0; } // CurrentSource1 has 1.0 A
+ }
+ }
+ { Name Voltage_Cir ;
+ Case {
+ { Region VoltageSource1; Value 10.0; } // VoltageSource1 has 10.0 V
+ }
+ }
+ { Name ElectricalCircuit ; Type Network ;
+ Case Circuit1 {
+ // Here the circuit netlist is implemented. The important thing is that
+ // each surface of the 3D geometry which is connected to the circuit
+ // (Sur_Ele) has a single voltage and current, but in the netlist they
+ // appear with two terminals (two nodes). The voltage between the two
+ // nodes corresponds to the absolute voltage of the surface. The current
+ // through the two nodes is the current flowing through the surface. In
+ // this simple case here we need the voltage drop across a cube (e.g. cube
+ // 1) in the netlist rather than the absolute voltage. Therefore the
+ // bottom-voltage has to be subtracted from the top-voltage. This is
+ // accomplished by connecting both surfaces anti-serial (e.g. at node 4;
+ // The voltage drop across cube 1 is then between node 2 and node 10).
+
+ { Region VoltageSource1; Branch{1, 10}; }
+ { Region Res1; Branch{1, 2}; }
+ // Voltage between node 2 and node 4 corresponds to the absolute voltage
+ // of Top1:
+ { Region Top1; Branch{2, 4}; }
+ // Note the reverse order of the nodes here - it's for subtracting the
+ // bottom voltage:
+ { Region Bottom1; Branch{10, 4}; }
+ // Voltage between node 3 and node 5 corresponds to the absolute voltage
+ // of Top2:
+ { Region Res2; Branch{2, 3}; }
+ { Region Top2; Branch{3, 5}; }
+ // Note the reverse order of the nodes here! It's for subtracting the
+ // bottom voltage.:
+ { Region Bottom2; Branch{10, 5}; }
+
+ { Region CurrentSource1; Branch{6, 10}; }
+ { Region Res3; Branch{6, 10}; }
+ // Voltage between node 6 and node 8 corresponds to the absolute voltage
+ // of Top3:
+ { Region Top3; Branch{6, 8}; }
+ // Note the reverse order of the nodes here! It's for subtracting the
+ // bottom voltage.:
+ { Region Bottom3; Branch{10, 8}; }
+ { Region Res4; Branch{6, 7}; }
+ // Voltage between node 7 and node 9 corresponds to the absolute voltage
+ // of Top4:
+ { Region Top4; Branch{7, 9}; }
+ // Note the reverse order of the nodes here! It's for subtracting the
+ // bottom voltage:
+ { Region Bottom4; Branch{10, 9}; }
+ }
+ }
+}
+
+Jacobian {
+ { Name JVol ;
+ Case {
+ { Region Region[{Vol_Ele}]; Jacobian Vol; }
+ { Region Region[{Sur_Ele}]; Jacobian Sur; }
+ }
+ }
+}
+
+Integration {
+ { Name I1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 5 ; }
+ { GeoElement Triangle ; NumberOfPoints 6 ; }
+ { GeoElement Quadrangle ; NumberOfPoints 7 ; }
+ { GeoElement Tetrahedron ; NumberOfPoints 15 ; }
+ { GeoElement Hexahedron ; NumberOfPoints 34 ; }
+ { GeoElement Prism ; NumberOfPoints 21 ; }
+ }
+ }
+ }
+ }
+}
+
+
+FunctionSpace {
+ { Name Hgrad_V; Type Form0;
+ BasisFunction {
+ // All nodes but the grouped nodes of the surfaces for connecting the
+ // circuitry
+ { Name sn; NameOfCoef vn; Function BF_Node;
+ Support VolWithSur_Ele; Entity NodesOf[ All, Not Sur_Ele ]; }
+ // Grouped nodes: Each surface of Sur_Ele has just one single voltage ->
+ // each surface is just one node effectively
+ { Name sf; NameOfCoef vf; Function BF_GroupOfNodes;
+ Support VolWithSur_Ele; Entity GroupsOfNodesOf[ Sur_Ele ]; }
+ }
+ GlobalQuantity {
+ // Voltage of the surfaces of Sur_Ele
+ { Name U; Type AliasOf; NameOfCoef vf; }
+ // Current flowing through the surfaces
+ { Name I; Type AssociatedWith; NameOfCoef vf; }
+ }
+ Constraint {
+ { NameOfCoef vn; EntityType NodesOf ; NameOfConstraint Voltage_3D; }
+ { NameOfCoef U; EntityType GroupsOfNodesOf ; NameOfConstraint Voltage_3D; }
+ { NameOfCoef I; EntityType GroupsOfNodesOf ; NameOfConstraint Current_3D; }
+ }
+ }
+
+ { Name Hregion_Cir; Type Scalar;
+ BasisFunction {
+ { Name sn; NameOfCoef ir; Function BF_Region;
+ Support Domain_Cir; Entity Domain_Cir; }
+ }
+ GlobalQuantity {
+ { Name Iz; Type AliasOf; NameOfCoef ir; }
+ { Name Uz; Type AssociatedWith; NameOfCoef ir; }
+ }
+ Constraint {
+ { NameOfCoef Uz ; EntityType Region ; NameOfConstraint Voltage_Cir ; }
+ { NameOfCoef Iz ; EntityType Region ; NameOfConstraint Current_Cir ; }
+ }
+ }
+}
+
+Formulation {
+ { Name Jstatic; Type FemEquation;
+ Quantity {
+ { Name V; Type Local; NameOfSpace Hgrad_V; }
+ { Name U; Type Global; NameOfSpace Hgrad_V [U]; }
+ { Name I; Type Global; NameOfSpace Hgrad_V [I]; }
+ { Name Uz; Type Global; NameOfSpace Hregion_Cir [Uz]; }
+ { Name Iz; Type Global; NameOfSpace Hregion_Cir [Iz]; }
+ }
+ Equation {
+ // FEM domain
+ Galerkin { [ kappa[] * Dof{d V}, {d V} ];
+ In Vol_Ele; Integration I1; Jacobian JVol; }
+ GlobalTerm { [ Dof{I} , {U} ]; In Sur_Ele; }
+
+ // Circuit related terms: Resistance equation
+ GlobalTerm { [ Dof{Uz}, {Iz} ]; In Resistance_Cir; }
+ GlobalTerm { [ R[] * Dof{Iz}, {Iz} ]; In Resistance_Cir; }
+
+ // Here is the connection between the circuit and the finite element simulation
+ GlobalEquation { Type Network; NameOfConstraint ElectricalCircuit;
+ { Node {I}; Loop {U}; Equation {I}; In Sur_Ele; }
+ { Node {Iz}; Loop {Uz}; Equation {Uz}; In Domain_Cir; }
+ }
+ }
+ }
+}
+
+Resolution {
+ { Name Jstatic;
+ System {
+ { Name Jstatic; NameOfFormulation Jstatic; }
+ }
+ Operation {
+ InitSolution[Jstatic];
+ Generate[Jstatic];
+ Solve[Jstatic];
+ SaveSolution[Jstatic];
+ }
+ }
+}
+
+PostProcessing {
+ { Name Ele; NameOfFormulation Jstatic;
+ Quantity {
+ { Name V; Value{ Local{ [ {V} ]; In Vol_Ele; Jacobian JVol;} } }
+ { Name Jv; Value{ Local{ [ -kappa[]*{d V} ]; In Vol_Ele; Jacobian JVol;} } }
+ { Name J; Value{ Local{ [ Norm[kappa[]*{d V}] ]; In Vol_Ele; Jacobian JVol;} } }
+ { Name U; Value { Term { [ {U} ]; In Sur_Ele; } } }
+ { Name I; Value { Term { [ {I} ]; In Sur_Ele; } } }
+ }
+ }
+}
+
+PostOperation {
+ // Absolute voltage [V]
+ { Name V; NameOfPostProcessing Ele;
+ Operation {
+ Print[ V, OnElementsOf Vol_Ele, File "Result_V.pos"];
+ }
+ }
+ // Current density vectors [A/m^2]
+ { Name J_vectors; NameOfPostProcessing Ele;
+ Operation {
+ Print[ Jv, OnElementsOf Vol_Ele, File "Result_J_vectors.pos"];
+ }
+ }
+ // Magnitude of current density [A/m^2]
+ { Name J_magnitude; NameOfPostProcessing Ele;
+ Operation {
+ Print[ J, OnElementsOf Vol_Ele, File "Result_J.pos"];
+ }
+ }
+ // Store the results in an ASCII file
+ { Name Cube_Top_Values_as_ASCII_File; NameOfPostProcessing Ele;
+ Operation {
+ Print[U, OnRegion Top1, File "Result_CubeTopValues.txt", Format SimpleTable];
+ Print[U, OnRegion Top2, File > "Result_CubeTopValues.txt", Format SimpleTable];
+ Print[U, OnRegion Top3, File > "Result_CubeTopValues.txt", Format SimpleTable];
+ Print[U, OnRegion Top4, File > "Result_CubeTopValues.txt", Format SimpleTable];
+
+ Print[I, OnRegion Top1, File > "Result_CubeTopValues.txt", Format SimpleTable];
+ Print[I, OnRegion Top2, File > "Result_CubeTopValues.txt", Format SimpleTable];
+ Print[I, OnRegion Top3, File > "Result_CubeTopValues.txt", Format SimpleTable];
+ Print[I, OnRegion Top4, File > "Result_CubeTopValues.txt", Format SimpleTable];
+ }
+ }
+}
diff --git a/CircuitCoupling/R_circuit_common.pro b/CircuitCoupling/R_circuit_common.pro
new file mode 100644
index 0000000000000000000000000000000000000000..bb6a5d3bac6def106a8ba8714293dd2ac8927171
--- /dev/null
+++ b/CircuitCoupling/R_circuit_common.pro
@@ -0,0 +1,15 @@
+// Physical numbers for volumes
+Cube1 = 1001;
+Cube2 = 1002;
+Cube3 = 1003;
+Cube4 = 1004;
+
+// Physical numbers for surfaces
+Cube1Top = 2001;
+Cube1Bottom = 2002;
+Cube2Top = 2003;
+Cube2Bottom = 2004;
+Cube3Top = 2005;
+Cube3Bottom = 2006;
+Cube4Top = 2007;
+Cube4Bottom = 2008;
\ No newline at end of file