Apparent discrepancy between website instructions and performance of "Specify Coefficients" approaches

I am working with a PDE that has one independent variable and varies only in space(m = 0, d = 0, only one equation). I have a 2D domain with multiple subdomains. Every subdomain has its own c, a, and f coefficients. At first, I used the specifyCoefficients() function in a loop over ever face/subdomain, assigning to each one its own constant coefficients. Now I am attempting to use a function to define those same coefficients, but the number of regions reported by the coefficients (using disp(length(region.x)) inside each function) is not consistent between c, f, and a, and is not consistent with the number suggested by the MathWorks website (see below). The contour plots of the two different solutions are attached (the only difference between them is the method I used to assign "a"), and below I will list the specifics of my pde object, as well as relevant code snippets.

Model specifics:

Number of faces = 200

Number of mesh nodes = 4890

Number of elements = 9532 ('linear' mesh)

Code used when assigning coefficients:

% Geometry and Mesh Specification
[geom_setup,geom_bt] = decsg(gdm,g_formula,g_name);
geom_return = geometryFromEdges(Model,geom_setup);%updates PDE model with geometry
generateMesh(Model,'GeometricOrder','linear');
[point_info,edge_info,tri_info] = meshToPet(Model.Mesh);
f_sub_num = tri_info(end,:)';%save face numbers associated with mesh triangles
f_elements = length(f_sub_num);%save number of triangles

% MAIN LOOP 1 - Assign coefficients for convection 'f'
for i = 1:f_elements

      % Obtain face number for mesh element
      face_num = f_sub_num(i);

      % Obtain voxel number given face number
      voxel_num = find(geom_bt(face_num,:),1);%this comes from decsg

      % Obtain location in data matrix from voxel number
      ii = g_pts_y - (mod(voxel_num,g_pts_y) +...
          g_pts_y*(mod(voxel_num,g_pts_y) == 0)) + 1;
      jj = floor((voxel_num-1)/g_pts_y) + 1;

      %Corresponding column of CM
      CM_col = CM(:,data(ii,jj)==CM(1,:));%matrix that I used as a look-up table for coefficient values

      %Read the values
      f_sub_coef(i) = CM_col(4);%global variable called inside @F_coefficient
  end

% MAIN LOOP 2 - Assign coefficients for absorption 'a' and diffusion 'c'
% Exactly 200 items are generated within Model.EquationCoefficients.CoefficientAssignments as expected
for i = 1:g_pts_y
    for j = 1:g_pts_x

          %Determine which voxel to assign coefficient to
          voxel_num = (g_pts_y-(i-1)) + g_pts_y*(j-1);

          %Corresponding column of CM
          CM_col = CM(:,data(i,j)==CM(1,:));%look up the correct values in table

          %Read the values
          c = CM_col(2);
          a = CM_col(3);
          f = @F_coeffunction;

          %Determine the corresponding face number
          face_num = find(geom_bt(:,voxel_num),1);

          %Assign coefficient value
          specifyCoefficients(Model,'m',m,'d',d,...
              'c',c,'a',a,'f',f,'Face',face_num);
      end
  end

%SEPARATE FILE - f function - simplified for discussion
function [f] = F_coeffunction(region,state)

global f_sub_coef
disp(['F length ' num2str(length(region.subdomain))]);
f = f_sub_coef(1:length(region.subdomain));

end

Note that length(region.x) = 9532 (the number of elements). However, if I attempt to use the same code seen in F_coefficient, for the "a" coefficient (the only difference is the value read from the CM_col plus an extra global variable so that "f" and "a" do not overwrite each other's values), I run into the following problem. When length(region.x) is called inside of a, it returns 3x9532, which corresponds to numel(Model.Mesh.Elements), not length(Model.Mesh.Elements(1,:))! The website's instructions say "Nr is equal to the length of the region.x or any other region field" for ALL coefficients a, c, and f.

Also, while running solvepde(Model), I noticed that length(region.x) changes from 1 to 9532. I thought that region.x was strictly related to the number of elements in the mesh, which caused problems when I tried to assign coefficients to each triangle as shown above. This led to the patching the code using f = "f_sub_coef(1:length(region.subdomain));".

QUESTION 1: Why does region.x change size for the f coefficient during the execution of solvepde(Model), and given that variation, how do I assign f coefficients for each triangle?

QUESTION 2: Why does region.x change size for "a" and "f"?

Secondly, my goal is to assign different a, c, and f values for each subdomain. It is easy to go back and forth between mesh triangles and face numbers using the outputs of meshToPet(). However, I have tried creating a coefficient vector for a using the code,

a_sub_coef = zeros(3,f_elements);%located above MAIN LOOP 1 -> global variable used in @A_coeff... as in @F_coeff...
a_sub_coef = reshape(a_sub_coef,1,numel(a_sub_coef));%located below MAIN LOOP 1

This way, there is a 1-1 correspondence between the columns of a_sub_coef and the triangles in the Mesh. I then assume that, after reusing reshape(), there is a 1-1 correspondence between the elements of a_sub_coef and region.x. However, when executing this code I obtain results that are different from the first approach.

QUESTION 3: How should I write @A_coefficient(region,state) so that I perform an assignment in the exact same fashion as looping over the faces and using specifyCoefficients(...)?

Recall, each face (therefore, the mesh triangles corresponding to those faces) should each have their own constant values of a, c, and f.

Thank you for your help with this very involved question!!

Hyperlinks in question: https://www.mathworks.com/help/pde/ug/specifycoefficients.html

https://www.mathworks.com/help/pde/ug/c-coefficient-for-systems-for-specifycoefficients.html

https://www.mathworks.com/help/pde/ug/f-coefficient-for-specifycoefficients.html

https://www.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html

Answer 1

[Steven Lord]

The online documentation is for the most recent release (currently release R2017a.) Are you using an older release? If so what happens when you use the versions of those pages included in the locally installed documentation for the release you're using?

[Eric]

Thank you for your quick response. I am using R2017a. The documentation with my local copy appears to be the same as the documentation online.