First version of systematic bases library for the hexahedron submitte…

…d to TAP

+assembly/getGlobalMatrices.m

@@ -0,0 +1,58 @@
+function [FEMatrices] = getGlobalMatrices(meshObj)
+    % GETGLOBALMATRICES computes the stiffness and mass matrices for the
+    % generalized eigenvalue problem.
+    % GETGLOBALMATRICES(meshObj) takes the mesh passed as parameter and
+    % returns a structure with two fields: stiffness and mass.
+
+    % Preallocation of sparse matrices.
+    nonZeros = meshObj.getNNZ();
+    irnValues = zeros(nonZeros,1);
+    jcnValues = zeros(nonZeros,1);
+    aValues = zeros(nonZeros,1);
+    a2Values = zeros(nonZeros,1);
+    sparseIndex = 0;
+
+
+    % We load here the basis functions computed before.
+    refBasesObj = bases.SystematicBases();
+    percentage = round(length(meshObj.elements)/10);
+
+    for elemIndex = 1:length(meshObj.elements)
+
+        % We obtain the different objects that represent the object and the
+        % material.
+        currentElemObj = meshObj.elements(elemIndex);
+        solidObj = meshObj.solids(currentElemObj.solidId);
+        % We compute the straight Jacobian Matrix.
+        jacobianMatrix = currentElemObj.computeJacobianMatrix();
+
+        % We compute the stiffness and mass matrices per element.
+        stiffElementMatrix = assembly.getStiffMatrix(refBasesObj, currentElemObj, jacobianMatrix, solidObj.materialObj.muRelative);
+        massElementMatrix = assembly.getMassMatrix(refBasesObj, currentElemObj, jacobianMatrix, solidObj.materialObj.epsilonRelative);
+
+        % We get the whole array of dofs to perform the finite element assembly.
+        dofIndices = currentElemObj.getDOF();
+
+        % It is stored columnwise (irn are rows, jcn are columns).
+        irnValues (sparseIndex+1:sparseIndex+length(dofIndices)^2) = reshape(dofIndices*ones(1,length(dofIndices)),[],1);
+        jcnValues (sparseIndex+1:sparseIndex+length(dofIndices)^2) = reshape(ones(length(dofIndices),1)*dofIndices',[],1);
+
+        % We impose the dirichlet boundary conditions to each element matrix.
+        localDirichletMatrix = assembly.imposeLocalDirichletMatrix (meshObj, currentElemObj, stiffElementMatrix);
+        aValues(sparseIndex+1:sparseIndex+length(dofIndices)^2) = reshape(localDirichletMatrix,[],1);
+
+        localDirichletMatrix = assembly.imposeLocalDirichletMatrix (meshObj, currentElemObj, massElementMatrix);
+        a2Values(sparseIndex+1:sparseIndex+length(dofIndices)^2) = reshape(localDirichletMatrix,[],1);
+
+        % To help with sparse indexing.
+        sparseIndex = sparseIndex + length(dofIndices)^2;
+
+        % For displaying purposes.
+        if (elemIndex/percentage == round(elemIndex/percentage))
+            disp([num2str(round(elemIndex/percentage)*10),  '% assembly done'])
+        end
+    end
+
+    % We perform here the finite element assembly.
+    FEMatrices.stiffness = sparse(irnValues(1:sparseIndex), jcnValues(1:sparseIndex), aValues(1:sparseIndex));
+    FEMatrices.mass = sparse(irnValues(1:sparseIndex), jcnValues(1:sparseIndex), a2Values(1:sparseIndex));

+assembly/getMassMatrix.m

@@ -0,0 +1,33 @@
+function massMatrix = getMassMatrix(refBasesObj, elemObj, jacobianMatrix, materialPropertyObj)
+    % GETMASSMATRIX - returns the mass matrix implemented with numerical
+    % integration.
+    % GETMASSMATRIX(refBasesObj, elemObj, jacobianMatrix, materialPropertyObj)
+    % takes the basis functions to use (defined in the reference element), the
+    % element where we perform the numerical integration for the mass matrix,
+    % the jacobian matrix to map the basis function to the real element,
+    % and the material to get the inverse of the magnetic permeability.
+
+    % We get the gauss points to perform the numerical integration.
+    integrationOrder = 2*elemObj.pOrder;
+    [gaussWeights,gaussPoints] = math.getGaussPoints(integrationOrder);
+
+    % Initialization.
+    massMatrix = zeros(54);
+
+    % This is assumed to be constant within the element (here we don't use
+    % curved or inhomogeneous materials).
+    epsilonR = materialPropertyObj.evaluateProperty();
+    detJacob = det(jacobianMatrix);
+
+    % We need to run through all the points for the numerical integration.
+    for gaussIndex = 1:length(gaussWeights)
+
+        % We evaluate the curl of the polynomials in the reference element, and
+        % then map it to the real element.
+        allBasesEvaluated = refBasesObj.getBasesEvaluated(gaussPoints(gaussIndex,:));
+        allBasesReal = jacobianMatrix\allBasesEvaluated.';
+
+        % We accumulate to perform the numerical integration.
+        massMatrix = massMatrix + gaussWeights(gaussIndex)*allBasesReal.'*epsilonR*(allBasesReal)*abs(detJacob);
+
+    end

+assembly/getStiffMatrix.m

@@ -0,0 +1,32 @@
+function stiffMatrix = getStiffMatrix(refBasesObj, elemObj, jacobianMatrix, materialPropertyObj)
+    % GETSTIFFMATRIX - returns the stiffness matrix implemented with numerical
+    % integration.
+    % GETSTIFFMATRIX(refBasesObj, elemObj, jacobianMatrix, materialPropertyObj)
+    % takes the basis functions to use (defined in the reference element), the
+    % element where we perform the numerical integration for the stiffness matrix,
+    % the jacobian matrix to map the basis function to the real element,
+    % and the material to get the inverse of the magnetic permeability.
+
+    % We get the gauss points to perform the numerical integration.
+    integrationOrder = 2*elemObj.pOrder;
+    [gaussWeights,gaussPoints] = math.getGaussPoints(integrationOrder);
+
+    % Initialization.
+    stiffMatrix = zeros(54);
+
+    % This is assumed to be constant within the element (here we don't use
+    % curved or inhomogeneous materials).
+    invMuRelative = materialPropertyObj.evaluateProperty();
+    detJacob = det(jacobianMatrix);
+
+    % We need to run through all the points for the numerical integration.
+    for gaussIndex = 1:length(gaussWeights)
+
+        % We evaluate the curl of the polynomials in the reference element, and
+        % then map it to the real element.
+        allCurlBasesEvaluated = refBasesObj.getCurlBasesEvaluated(gaussPoints(gaussIndex,:));
+        allCurlBasesReal = allCurlBasesEvaluated*jacobianMatrix/detJacob;
+
+        % We accumulate to perform the numerical integration.
+        stiffMatrix = stiffMatrix + gaussWeights(gaussIndex)*allCurlBasesReal*invMuRelative*(allCurlBasesReal.')*abs(detJacob);
+    end

+assembly/imposeLocalDirichletMatrix.m

@@ -0,0 +1,53 @@
+function [dirichletMatrix] = imposeLocalDirichletMatrix(meshObj, elemObj, matrix)
+    % IMPOSELOCALDIRICHLET returns an element matrix where we have set to 1 the diagonal and
+    % removed the remaining terms in row and columns.
+    % IMPOSELOCALDIRICHLET(meshObj, elemObj, matrix) takes the boundary conditions from the
+    % mesh passed as parameter, runs through all the edges and faces from elemObj and
+    % apply the Dirichlet boundary conditions (setting to 0 the rows and columns and to
+    % one the diagonals where Dirichlet is applied) to the matrix passed as parameter.
+
+    % Counter to check all the dofs.
+    currentDOF = 0;
+    dirichletMatrix = matrix;
+
+    % We run through all the edges and check if the edge belongs to a Dirichlet
+    % boundary condition.
+    for indexEdge = 1:length(elemObj.edges)
+        polyFaceId = elemObj.edges(indexEdge).polyFaceId;
+        lengthDOF = length(elemObj.edges(indexEdge).dof);
+
+        [dirichletMatrix] = checkDirichlet (dirichletMatrix, meshObj, currentDOF, lengthDOF, polyFaceId);
+        currentDOF = currentDOF + lengthDOF;
+    end
+
+    % The same for the faces.
+    for indexFace = 1:length(elemObj.faces)
+        polyFaceId = elemObj.faces(indexFace).polyFaceId;
+        lengthDOF = length(elemObj.faces(indexFace).dof);
+
+        [dirichletMatrix] = checkDirichlet (dirichletMatrix, meshObj, currentDOF, lengthDOF, polyFaceId);
+        currentDOF = currentDOF + lengthDOF;
+    end
+
+end
+
+function [dirichletMatrix] = checkDirichlet (dirichletMatrix, meshObj, currentDOF, lengthDOF, polyFaceId)
+    % Internal function to check for both entities (edge and face) if we need to apply Dirichlet.
+
+    dirichletString = "PerfectE";
+    if (polyFaceId ~= 0)
+        if (strcmp(meshObj.polyFaces(polyFaceId).polyFaceType, dirichletString))
+            dirichletMatrix = imposeDirichletInRowAndCol(dirichletMatrix, currentDOF, lengthDOF);
+        end
+    end
+
+end
+
+function [dirichletMatrix] = imposeDirichletInRowAndCol (dirichletMatrix, currentDOF, lengthDOF)
+    % Internal function to remove the Dirichlet unknown.
+
+    dirichletMatrix(currentDOF+1:currentDOF+lengthDOF,:) = 0;
+    dirichletMatrix(:,currentDOF+1:currentDOF+lengthDOF) = 0;
+    dirichletMatrix(currentDOF+1:currentDOF+lengthDOF,currentDOF+1:currentDOF+lengthDOF) = eye(lengthDOF);
+
+end

+bases/@SystematicBases/SystematicBases.m

@@ -0,0 +1,184 @@
+classdef SystematicBases < handle
+
+    % Class that contains the mixed-order curl-conforming basis functions
+    % presented in the paper.
+
+    properties
+
+        % Reference elements where the functions belong.
+        refElement geometry.Element
+
+        % Order of the basis functions (only p=2 supported).
+        pOrder double
+
+        % Array of all the polynomials that compose the mixed-order.
+        % space of functions.
+        polyAllFunctions(:,3) math.Polynomial
+
+        % Array of the curl of polynomials in polyAllFunctions.
+        polyAllCurlFunctions(:,3) math.Polynomial
+
+        % Symbolic representation of the polynomials.
+        symAllFunctions(:,3) sym
+
+        % Symbolic representation of the curl of the polynomials.
+        symAllCurlFunctions(:,3) sym
+
+        % Array of symbolic coordinates.
+        realCoordinates math.RealCoordinates
+
+        % Reference q polynomials for the edge dof (6).
+        qRefEdge math.Polynomial
+        % Reference q polynomials for the face dof (7).
+        qRefRectangular (:,:,3) math.Polynomial
+        % Reference q polynomials for the volume dof (8).
+        qRefHexa (:,3) math.Polynomial
+
+        % Reference coefficients obtained when executing getCoefficients().
+        referenceCoefficients(:,:) double
+
+    end
+
+    methods
+        function systematicBasesObject = SystematicBases()
+            % Constructor of the basis functions. This method gets the
+            % saved object from precomputedBases, or initializes all the
+            % machinery to obtain the coefficients solving the dual bases.
+
+            basesName = 'systematicBasesHexa2';
+            % Get the precomputed bases from memory.
+            if (isfile(strcat('+bases/precomputedBases/',basesName,'.mat')))
+                load(strcat('+bases/precomputedBases/',basesName,'.mat'), 'tempBases')
+
+                systematicBasesObject.refElement = tempBases.refElement;
+                systematicBasesObject.pOrder = tempBases.pOrder;
+                systematicBasesObject.polyAllFunctions = tempBases.polyAllFunctions;
+                systematicBasesObject.polyAllCurlFunctions = tempBases.polyAllCurlFunctions;
+                systematicBasesObject.symAllFunctions = tempBases.symAllFunctions;
+                systematicBasesObject.symAllCurlFunctions = tempBases.symAllCurlFunctions;
+                systematicBasesObject.realCoordinates = tempBases.realCoordinates;
+                systematicBasesObject.qRefEdge = tempBases.qRefEdge;
+                systematicBasesObject.qRefRectangular = tempBases.qRefRectangular;
+                systematicBasesObject.qRefHexa = tempBases.qRefHexa;
+                systematicBasesObject.referenceCoefficients = tempBases.referenceCoefficients;
+            else
+                % We obtain the reference element defined in Fig. 1.
+                systematicBasesObject.refElement = geometry.Element ('referenceElement');
+                % We define the order, which is important for assigning the dofs.
+                systematicBasesObject.pOrder = 2;
+                % We obtain symbolic variables [x, y, z].
+                systematicBasesObject.realCoordinates = math.RealCoordinates(true);
+                xyz = systematicBasesObject.realCoordinates.coordinates;
+
+                % These are the symbolic representation of the polynomials that
+                % compose the space of functions.
+                symAllFunctions = bases.getPolynomials(systematicBasesObject.pOrder);
+                lengthPolymials = size(symAllFunctions,1);
+
+                % We obtain the curl of all these polynomials.
+                symAllCurlFunctions = sym(zeros(size(symAllFunctions)));
+                for indexPolymial = 1:lengthPolymials
+                    symAllCurlFunctions(indexPolymial,:) = curl(symAllFunctions(indexPolymial,:), xyz).';
+                end
+                systematicBasesObject.symAllFunctions = symAllFunctions;
+                systematicBasesObject.symAllCurlFunctions = symAllCurlFunctions;
+
+                % Now, we get the polynomial representation of these polynomial.s
+                systematicBasesObject.polyAllFunctions(lengthPolymials, 3) = math.Polynomial;
+                systematicBasesObject.polyAllCurlFunctions(lengthPolymials, 3) = math.Polynomial;
+                for indexPolymial = 1:lengthPolymials
+                    for indexDimension = 1:3
+                        systematicBasesObject.polyAllFunctions(indexPolymial, indexDimension) = math.Polynomial(systematicBasesObject.pOrder*ones(1,3), symAllFunctions(indexPolymial, indexDimension));
+                        systematicBasesObject.polyAllCurlFunctions(indexPolymial, indexDimension) = math.Polynomial(systematicBasesObject.pOrder*ones(1,3), symAllCurlFunctions(indexPolymial, indexDimension));
+                    end
+                end
+
+                % To save space in the following.
+                t1 = systematicBasesObject.realCoordinates.coordinates(1);
+                t2 = systematicBasesObject.realCoordinates.coordinates(2);
+                t3 = systematicBasesObject.realCoordinates.coordinates(3);
+
+                % We define the same q for all the edges to discretize (6).
+                % All these q are contained in Table I.
+                qRefEdge(2) = math.Polynomial();
+                qRefEdge(1) = math.Polynomial(systematicBasesObject.pOrder*ones(1,3), 1-t1);
+                qRefEdge(2) = math.Polynomial(systematicBasesObject.pOrder*ones(1,3), t1);
+
+                % The vector polynomials in (7) are defined here for all the
+                % faces following the order in the reference element given in
+                % getFaceByLocalEdges.
+                qSymRectangular = sym(zeros(6,4,3));
+                qSymRectangular(1,1,:) = [1-t1 0 0];
+                qSymRectangular(1,2,:) = [t1 0 0];
+                qSymRectangular(1,3,:) = [0 1-t2 0];
+                qSymRectangular(1,4,:) = [0 t2 0];
+                qSymRectangular(2,1,:) = [1-t1 0 0];
+                qSymRectangular(2,2,:) = [t1 0 0];
+                qSymRectangular(2,3,:) = [0 1-t2 0];
+                qSymRectangular(2,4,:) = [0 t2 0];
+                qSymRectangular(3,1,:) = [1-t1 0 0];
+                qSymRectangular(3,2,:) = [t1 0 0];
+                qSymRectangular(3,3,:) = [0 0 1-t3];
+                qSymRectangular(3,4,:) = [0 0 t3];
+                qSymRectangular(4,1,:) = [0 1-t2 0];
+                qSymRectangular(4,2,:) = [0 t2 0];
+                qSymRectangular(4,3,:) = [0 0 1-t3];
+                qSymRectangular(4,4,:) = [0 0 t3];
+                qSymRectangular(5,1,:) = [1-t1 0 0];
+                qSymRectangular(5,2,:) = [t1 0 0];
+                qSymRectangular(5,3,:) = [0 0 1-t3];
+                qSymRectangular(5,4,:) = [0 0 t3];
+                qSymRectangular(6,1,:) = [0 1-t2 0];
+                qSymRectangular(6,2,:) = [0 t2 0];
+                qSymRectangular(6,3,:) = [0 0 1-t3];
+                qSymRectangular(6,4,:) = [0 0 t3];
+
+                % We obtain now the polynomial representation.
+                qRefRectangular(6,4,3) = math.Polynomial();
+                for indexFace = 1:6
+                    for indexComponent = 1:3
+                        for indexQ = 1:4
+                            qRefRectangular(indexFace,indexQ, indexComponent) = math.Polynomial(systematicBasesObject.pOrder*ones(1,3), qSymRectangular(indexFace, indexQ, indexComponent));
+                        end
+                    end
+                end
+
+                % Finally, we define the vector polynomial q needed in (8)
+                qSymHexa = sym(zeros(6,3));
+                qSymHexa(1,:) = [1-t1 0 0];
+                qSymHexa(2,:) = [t1 0 0];
+                qSymHexa(3,:) = [0 1-t2 0];
+                qSymHexa(4,:) = [0 t2 0];
+                qSymHexa(5,:) = [0 0 1-t3];
+                qSymHexa(6,:) = [0 0 t3];
+
+                % Obtaining now the polynomial representation.
+                qRefHexa(6,3) = math.Polynomial();
+                for indexComponent = 1:3
+                    for indexQ = 1:6
+                        qRefHexa(indexQ, indexComponent) = math.Polynomial(systematicBasesObject.pOrder*ones(1,3), qSymHexa(indexQ, indexComponent));
+                    end
+                end
+
+                systematicBasesObject.qRefEdge = qRefEdge;
+                systematicBasesObject.qRefRectangular = qRefRectangular;
+                systematicBasesObject.qRefHexa = qRefHexa;
+            end
+
+        end
+    end
+
+    methods (Static)
+        function removeBases()
+            % This method removes the copy stored in memory to generate a new
+            % set of basis functions.
+
+            if (isfile('+bases/precomputedBases/systematicBasesHexa2.mat'))
+                delete '+bases/precomputedBases/systematicBasesHexa2.mat'
+            end
+
+        end
+
+    end
+
+end