diff --git a/smart/mesh_tools.py b/smart/mesh_tools.py
index 576a1be..e1a700a 100644
--- a/smart/mesh_tools.py
+++ b/smart/mesh_tools.py
@@ -38,7 +38,7 @@
 ]
 
 
-def implicit_curve(boundExpr):
+def implicit_curve(boundExpr, num_points=51):
     """
     Output (r,z) coordinates given a string expression in the r-z plane
     Currently the function only outputs cell coordinates in the right half
@@ -47,6 +47,8 @@ def implicit_curve(boundExpr):
     which corresponds to a substrate in most cases.
     The (r,z) coordinates start at the maximum z value for r = 0
     and end when either z = 0 or r = 0.
+    num_points can be increased for more accuracy in generating the contour,
+    but this function will run more slowly for very high num_points.
 
     Examples:
     * boundExpr = "r**2 + z**2 - 1": defines a quarter circle from (0, 1) to (1, 0)
@@ -57,13 +59,16 @@ def implicit_curve(boundExpr):
     z = sym.Symbol("z", real=True)
     outerExpr0 = outerExprRef.subs({"r": 0, "z": z})
     z0 = solveset_real(outerExpr0, z)
+    z0Array = np.array([float(z0Val) for z0Val in z0])
     rVals = [0.0]
     zVals = [float(max(z0))]
     rMax = solveset_real(outerExprRef.subs({"r": r, "z": 0.0}), r)
     if len(rMax) > 0:
-        sGap = max([float(max(z0)), float(max(rMax))]) / 50
+        sGap = max([float(max(z0)), float(max(rMax))]) / (num_points - 1)
+    elif len(z0) == 2 and np.all(z0Array > 0):  # then probably a nuclear contour
+        sGap = float(max(z0) - min(z0)) / (num_points - 1)
     else:
-        sGap = float(max(z0)) / 50
+        sGap = float(max(z0)) / (num_points - 1)
     curTan = [1, 0]
     while rVals[-1] >= 0 and zVals[-1] >= 0:
         rNext = rVals[-1] + curTan[0] * sGap
diff --git a/smart/model.py b/smart/model.py
index 809abf5..6e0cd21 100644
--- a/smart/model.py
+++ b/smart/model.py
@@ -632,7 +632,8 @@ def _init_2_4_check_for_unused_parameters_species_compartments(self):
                 raise ValueError(print_str)
 
     def _init_2_5_link_compartments_to_species(self):
-        """Linking compartments and compartment dimensionality to species"""
+        """Linking compartments and compartment dimensionality to species,
+        check for consistency of diffusion coefficient definition"""
         logger.debug(
             "Linking compartments and compartment dimensionality to species",
             extra=dict(format_type="log"),
@@ -640,6 +641,12 @@ def _init_2_5_link_compartments_to_species(self):
         for species in self.sc:
             species.compartment = self.cc[species.compartment_name]
             species.dimensionality = self.cc[species.compartment_name].dimensionality
+            # convert diffusion coeff to units consistent with mesh
+            diffusion_conversion = species.diffusion_units.to(
+                species.compartment.compartment_units**2 / unit.s
+            )
+            species.diffusion_units = species.compartment.compartment_units**2 / unit.s
+            species.D *= diffusion_conversion.magnitude
 
     def _init_2_6_link_species_to_compartments(self):
         """Links species to compartments - a species is considered to be
@@ -1123,11 +1130,6 @@ def _init_5_2_create_variational_forms(self):
                     extra=dict(format_type="log"),
                 )
             else:
-                if Dform_units != mass_form_units:  # unit conversion for consistency
-                    diffusion_conversion = species.diffusion_units.to(
-                        species.compartment.compartment_units**2 / unit.s
-                    )
-                    D *= diffusion_conversion.magnitude
                 D_constant = d.Constant(D, name=f"D_{species.name}")
                 if self.config.flags["axisymmetric_model"]:
                     Dform = x[0] * D_constant * d.inner(d.grad(u), d.grad(v)) * dx