diff --git a/+assembly/getGlobalMatrices.m b/+assembly/getGlobalMatrices.m new file mode 100644 index 0000000..cd3e856 --- /dev/null +++ b/+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)); diff --git a/+assembly/getMassMatrix.m b/+assembly/getMassMatrix.m new file mode 100644 index 0000000..0c722cf --- /dev/null +++ b/+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 \ No newline at end of file diff --git a/+assembly/getStiffMatrix.m b/+assembly/getStiffMatrix.m new file mode 100644 index 0000000..601eb07 --- /dev/null +++ b/+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 \ No newline at end of file diff --git a/+assembly/imposeLocalDirichletMatrix.m b/+assembly/imposeLocalDirichletMatrix.m new file mode 100644 index 0000000..23ce492 --- /dev/null +++ b/+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 \ No newline at end of file diff --git a/+bases/@SystematicBases/SystematicBases.m b/+bases/@SystematicBases/SystematicBases.m new file mode 100644 index 0000000..031e5d3 --- /dev/null +++ b/+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 \ No newline at end of file diff --git a/+bases/@SystematicBases/evaluateEdgeDOF.m b/+bases/@SystematicBases/evaluateEdgeDOF.m new file mode 100644 index 0000000..ef841b2 --- /dev/null +++ b/+bases/@SystematicBases/evaluateEdgeDOF.m @@ -0,0 +1,36 @@ +function edgeDOFs = evaluateEdgeDOF(sysBasesObj, edgeObj, functionToEvaluate) + % EVALUATEEDGEDOF returns the numerical integration of (6). + % EVALUATEEDGEDOF(sysBasesObj,edgeObj,functionToEvaluate) applys (6) at the + % edge of the element passed as parameter to the polynomial representation + % given in functionToEvaluate. This function returns the two dofs for the q + % polynomials defined in the SystematicBasis object. + + edgeUnitVector = edgeObj.getUnitVector(); + % We obtain the 1D jacobian for the edge. + [jacobianEdge, originPointInEdge] = edgeObj.getJacobian(); + % We need the length of the edge to evaluate the dof + lengthEdge = sqrt(det(jacobianEdge.'*jacobianEdge)); + + % We obtain the points for the numerical integration. + [gaussWeights, gaussPoints1D] = math.getGaussPointsForSegment(2*sysBasesObj.pOrder); + qArray = sysBasesObj.qRefEdge; + edgeDOFs = zeros(size(qArray)); + + % We apply the functional for all the q polynomials. + for indexQ = 1:length(qArray) + for indexPoint = 1:length(gaussWeights) + % To get the 3D representation of the gauss point for the segment. + point3D = jacobianEdge*gaussPoints1D(indexPoint)+originPointInEdge; + functionEvaluated = zeros(1,3); + % We evaluate the different polynomials. + for indexCoordinates = 1:3 + functionEvaluated(indexCoordinates) = functionToEvaluate(indexCoordinates).evaluatePolynomial(point3D); + end + qEvaluated = qArray(indexQ).evaluatePolynomial([gaussPoints1D(indexPoint) 0 0]); + % We perform here the numerical integration. + edgeDOFs(indexQ) = edgeDOFs(indexQ) + gaussWeights(indexPoint)*dot(functionEvaluated,edgeUnitVector)*qEvaluated*lengthEdge; + end + end + +end + diff --git a/+bases/@SystematicBases/evaluateFaceDOF.m b/+bases/@SystematicBases/evaluateFaceDOF.m new file mode 100644 index 0000000..3b8fbcb --- /dev/null +++ b/+bases/@SystematicBases/evaluateFaceDOF.m @@ -0,0 +1,36 @@ +function faceDOFs = evaluateFaceDOF(sysBasesObj, faceNumber, functionToEvaluate) + % EVALUATEFACEDOF returns the numerical integration of (7). + % EVALUATEFACEDOF(sysBasesOb,faceNumber,functionToEvaluate) applys + % (7) on the face of the element passed as parameter to the polynomial + % representation given in functionToEvaluate. This function returns the four + % dofs for the q polynomial vectors defined in the SystematicBasis object. + + % All the faces in the reference element have unit size. + detJacob2D = 1; + + % We obtain the gauss points to perform the numerical integration on a + % quadrilateral domain. + [gaussWeights, gaussPoints2D] = math.getGaussPointsForQuadrilateral(2*sysBasesObj.pOrder); + + qArray = squeeze(sysBasesObj.qRefRectangular(faceNumber,:,:)); + faceDOFs = zeros(size(qArray,1),1); + functionEvaluated = zeros(1,3); + qRefEvaluated = zeros(3,1); + + % We apply the functional for all the q polynomial vectors. + for indexQ = 1:size(qArray,1) + for indexPoint = 1:length(gaussWeights) + % To get the 3D representation of the gauss point for the segment. + gaussPoint3D = math.transform2Dto3DInReferenceElement(gaussPoints2D(indexPoint,:), faceNumber); + + % We evaluate the different polynomials. + for indexCoordinates = 1:3 + functionEvaluated(indexCoordinates) = functionToEvaluate(indexCoordinates).evaluatePolynomial(gaussPoint3D); + qRefEvaluated(indexCoordinates) = qArray(indexQ,indexCoordinates).evaluatePolynomial(gaussPoint3D); + end + % We perform here the numerical integration. + faceDOFs(indexQ) = faceDOFs(indexQ) + gaussWeights(indexPoint)*dot(functionEvaluated,qRefEvaluated)*detJacob2D; + end + end + +end diff --git a/+bases/@SystematicBases/evaluateVolumeDOF.m b/+bases/@SystematicBases/evaluateVolumeDOF.m new file mode 100644 index 0000000..d0d6a06 --- /dev/null +++ b/+bases/@SystematicBases/evaluateVolumeDOF.m @@ -0,0 +1,35 @@ +function volDOFs = evaluateVolumeDOF(sysBasesObj, elemObj, functionToEvaluate) + % EVALUATEVOLUMEDOF returns the numerical integration of (8). + % EVALUATEVOLUMEDOF(sysBasesObj,elemObj,functionToEvaluate) applys + % (8) for the interior of the element passed as parameter to the polynomial + % representation given in functionToEvaluate. This function returns the six + % dofs for the q polynomial vectors defined in the SystematicBasis object. + + % We obtain the volumetric jacobian for the element. + jacobianMatrix = elemObj.computeJacobianMatrix(); + + % We obtain the gauss points to perform the numerical integration on a + % hexahedral domain. + [gaussWeights,gaussPoints] = math.getGaussPoints(2*sysBasesObj.pOrder); + + qArray = sysBasesObj.qRefHexa; + volDOFs = zeros(size(qArray,1),1); + functionEvaluated = zeros(1,3); + qEvaluated = zeros(1,3); + + % We apply the functional for all the q polynomial vectors. + for indexQ = 1:size(qArray,1) + for indexPoint = 1:length(gaussWeights) + + % We evaluate all the polynomials involved. + for indexComponent = 1:3 + functionEvaluated(indexComponent) = functionToEvaluate(indexComponent).evaluatePolynomial(gaussPoints(indexPoint,:)); + % This should be used (probably) together with the unit vector + qEvaluated(indexComponent) = qArray(indexQ,indexComponent).evaluatePolynomial(gaussPoints(indexPoint,:)); + end + % We perform here the numerical integration. + volDOFs(indexQ) = volDOFs(indexQ) + gaussWeights(indexPoint)*dot(functionEvaluated,qEvaluated)*det(jacobianMatrix); + end + end + +end diff --git a/+bases/@SystematicBases/getBasesEvaluated.m b/+bases/@SystematicBases/getBasesEvaluated.m new file mode 100644 index 0000000..34a9fae --- /dev/null +++ b/+bases/@SystematicBases/getBasesEvaluated.m @@ -0,0 +1,25 @@ +function allBasesEvaluated = getBasesEvaluated(sysBasesObj, gaussPoint) + % GETBASESEVALUATED returns the curl of the 54 basis functions. + % GETBASESEVALUATED(sysBasesObj, gaussPoint) takes the linear combination of the + % curl of the polynomials stored in sysBasesObj and evaluate them in gaussPoint. Then, + % the linear combination with the coefficients obtained in GETCOEFFICIENTS is obtained. + % See also: GETCOEFFICIENTS + + coefficients = sysBasesObj.referenceCoefficients; + lengthBases = size(sysBasesObj.polyAllFunctions,1); + allBasesEvaluated = zeros(size(sysBasesObj.polyAllFunctions)); + polynomialsEvaluated = zeros(size(sysBasesObj.polyAllFunctions)); + + % We obtain the evaluation of the polynomial in the point passed as parameter. + for indexPolymials = 1:lengthBases + for indexCoordinates = 1:3 + polynomialsEvaluated(indexPolymials, indexCoordinates) = sysBasesObj.polyAllFunctions(indexPolymials, indexCoordinates).evaluatePolynomial(gaussPoint); + end + end + + % We apply the linear combination of the polynomials to obtain the 54 bases. + for indexBases = 1:lengthBases + allBasesEvaluated(indexBases,:) = coefficients(:,indexBases).'*polynomialsEvaluated; + end + +end \ No newline at end of file diff --git a/+bases/@SystematicBases/getCoefficients.m b/+bases/@SystematicBases/getCoefficients.m new file mode 100644 index 0000000..030345b --- /dev/null +++ b/+bases/@SystematicBases/getCoefficients.m @@ -0,0 +1,41 @@ +function coefficients = getCoefficients (sysBasesObj) + % GETCOEFFICIENTS Evaluates numerically the dofs in (6), (7), and (8) and + % obtains the coefficients of the polynomials contained in the mixed-order + % curl-conforming space through the dual basis imposition of these dofs. + % GETCOEFFICIENTS(sysBasesObj) takes an object of the SystematicBases class + % to provide the linear combination of the polynomials of the space. + % The result is shown in table II. + % See also GETPOLYNOMIALS, EVALUATEEDGEDOF, EVALUATEFACEDOF, EVALUATEVOLUMEDOF. + + % We use here always the reference element (it could also be done for the real element). + elemObj = sysBasesObj.refElement; + lengthBases = size(sysBasesObj.polyAllFunctions,1); + lhsMatrixGiMon = zeros(lengthBases); + + % For all the polynomials within the space... + for indexPolynomial = 1:lengthBases + % We get the polynomial representation (faster than symbolic variables). + polynomialToEvaluate = sysBasesObj.polyAllFunctions(indexPolynomial, :); + dofCounter = 0; + % We apply (6) at each edge. + for indexEdge = 1:length(elemObj.edges) + edgeObj = elemObj.edges(indexEdge); + edgeDOFs = sysBasesObj.evaluateEdgeDOF(edgeObj, polynomialToEvaluate); + lhsMatrixGiMon(dofCounter+1:dofCounter+length(edgeDOFs), indexPolynomial) = edgeDOFs; + dofCounter = dofCounter + length(edgeDOFs); + end + % We apply (7) on each face. + for indexFace = 1:length(elemObj.faces) + faceDOFs = sysBasesObj.evaluateFaceDOF(indexFace, polynomialToEvaluate); + lhsMatrixGiMon(dofCounter+1:dofCounter+length(faceDOFs), indexPolynomial) = faceDOFs; + dofCounter = dofCounter + length(faceDOFs); + end + % We apply (8) inside each element. + volumeDOFs = sysBasesObj.evaluateVolumeDOF(elemObj, polynomialToEvaluate); + lhsMatrixGiMon(dofCounter+1:dofCounter+length(volumeDOFs), indexPolynomial) = volumeDOFs; + end + + % We obtain the coefficients as the dual basis of the dofs. + coefficients = lhsMatrixGiMon\eye(lengthBases); + sysBasesObj.referenceCoefficients = coefficients; +end \ No newline at end of file diff --git a/+bases/@SystematicBases/getCurlBasesEvaluated.m b/+bases/@SystematicBases/getCurlBasesEvaluated.m new file mode 100644 index 0000000..f9b3fb5 --- /dev/null +++ b/+bases/@SystematicBases/getCurlBasesEvaluated.m @@ -0,0 +1,25 @@ +function allCurlBasesEvaluated = getCurlBasesEvaluated(sysBasesObj, gaussPoint) + % GETCURLBASESEVALUATED returns the curl of the 54 basis functions. + % GETCURLBASESEVALUATED(sysBasesObj, gaussPoint) takes the linear combination of the + % curl of the polynomials stored in sysBasesObj and evaluate them in gaussPoint. Then, + % the linear combination with the coefficients obtained in GETCOEFFICIENTS is obtained. + % See also: GETCOEFFICIENTS + + coefficients = sysBasesObj.referenceCoefficients; + lengthBases = size(sysBasesObj.polyAllFunctions,1); + allCurlBasesEvaluated = zeros(size(sysBasesObj.polyAllCurlFunctions)); + curlPolynomialsEvaluated = zeros(size(sysBasesObj.polyAllCurlFunctions)); + + % We obtain the curl of the polynomial in the point passed as parameter. + for indexPolymials = 1:lengthBases + for indexCoordinates = 1:3 + curlPolynomialsEvaluated(indexPolymials, indexCoordinates) = sysBasesObj.polyAllCurlFunctions(indexPolymials, indexCoordinates).evaluatePolynomial(gaussPoint); + end + end + + % We apply the linear combination of the polynomials to obtain the 54 bases. + for indexBases = 1:lengthBases + allCurlBasesEvaluated(indexBases,:) = coefficients(:,indexBases).'*curlPolynomialsEvaluated; + end + +end \ No newline at end of file diff --git a/+bases/@SystematicBases/saveCoefficients.m b/+bases/@SystematicBases/saveCoefficients.m new file mode 100644 index 0000000..498aa98 --- /dev/null +++ b/+bases/@SystematicBases/saveCoefficients.m @@ -0,0 +1,16 @@ +function saveCoefficients(systematicBasesObject) + % SAVECOEFFICIENTS store the whole object of systematic bases. + % SAVECOEFFICIENTS(systematicBasesObject) takes the object + % passed as parameter (including the coefficients obtained + % as dual bases) and store it to be faster in following calls. + + basesName = 'systematicBasesHexa2'; + tempBases = 'tempBases'; + functionsStruct.(tempBases) = systematicBasesObject; + fullFile = fullfile('+bases/precomputedBases/',basesName); + if (~exist('+bases/precomputedBases/', 'dir')) + mkdir ('+bases/precomputedBases/'); + end + save (fullFile, '-struct', 'functionsStruct'); + +end \ No newline at end of file diff --git a/+bases/@SystematicBases/showBases.m b/+bases/@SystematicBases/showBases.m new file mode 100644 index 0000000..b77bad1 --- /dev/null +++ b/+bases/@SystematicBases/showBases.m @@ -0,0 +1,98 @@ +function showBases(sysBasesObj, basesToShow) + % SHOWBASES displays the coefficients that multiply the polynomials + % contained in the spaace. + % SHOWBASES(sysBasesObj, basesToShow) uses the coefficients obtained + % from GETCOEFFICIENTS to display for the different bases the + % coefficients known as a1, a2... in (4). If basesToShow is omitted, + % all the functions are displayed. + + if (nargin == 1) + basesToShow = 1:54; + end + + tolerance = 1e-12; + + lengthBases = length(basesToShow); + + stringCoefficients = ''; + header = 'Ni '; + header2 = '----'; + for indexBases = 1:lengthBases + stringCoefficients = strcat(stringCoefficients,'|%4g'); + header = [header,'|%4i']; + header2 = [header2,'|----']; + end + stringCoefficients = strcat(stringCoefficients,'\n'); + header = [header,'\n']; + header2 = [header2,'\n']; + + nameCoefficients = {}; + nameCoefficients{1}= 'a1\t'; + nameCoefficients{2}= 'a2\t'; + nameCoefficients{3}= 'a3\t'; + nameCoefficients{4}= 'a4\t'; + nameCoefficients{5}= 'a5\t'; + nameCoefficients{6}= 'a6\t'; + nameCoefficients{7}= 'a7\t'; + nameCoefficients{8}= 'a8\t'; + nameCoefficients{9}= 'a9\t'; + nameCoefficients{10} = 'a10\t'; + nameCoefficients{11} = 'a11\t'; + nameCoefficients{12} = 'a12\t'; + nameCoefficients{13} = 'a13\t'; + nameCoefficients{14} = 'a14\t'; + nameCoefficients{15} = 'a15\t'; + nameCoefficients{16} = 'a16\t'; + nameCoefficients{17} = 'a17\t'; + nameCoefficients{18} = 'a18\t'; + nameCoefficients{19} = 'b1\t'; + nameCoefficients{20} = 'b2\t'; + nameCoefficients{21} = 'b3\t'; + nameCoefficients{22} = 'b4\t'; + nameCoefficients{23} = 'b5\t'; + nameCoefficients{24} = 'b6\t'; + nameCoefficients{25} = 'b7\t'; + nameCoefficients{26} = 'b8\t'; + nameCoefficients{27} = 'b9\t'; + nameCoefficients{28} = 'b10\t'; + nameCoefficients{29} = 'b11\t'; + nameCoefficients{30} = 'b12\t'; + nameCoefficients{31} = 'b13\t'; + nameCoefficients{32} = 'b14\t'; + nameCoefficients{33} = 'b15\t'; + nameCoefficients{34} = 'b16\t'; + nameCoefficients{35} = 'b17\t'; + nameCoefficients{36} = 'b18\t'; + nameCoefficients{37} = 'c1\t'; + nameCoefficients{38} = 'c2\t'; + nameCoefficients{39} = 'c3\t'; + nameCoefficients{40} = 'c4\t'; + nameCoefficients{41} = 'c5\t'; + nameCoefficients{42} = 'c6\t'; + nameCoefficients{43} = 'c7\t'; + nameCoefficients{44} = 'c8\t'; + nameCoefficients{45} = 'c9\t'; + nameCoefficients{46} = 'c10\t'; + nameCoefficients{47} = 'c11\t'; + nameCoefficients{48} = 'c12\t'; + nameCoefficients{49} = 'c13\t'; + nameCoefficients{50} = 'c14\t'; + nameCoefficients{51} = 'c15\t'; + nameCoefficients{52} = 'c16\t'; + nameCoefficients{53} = 'c17\t'; + nameCoefficients{54} = 'c18\t'; + + fprintf(header,basesToShow); + fprintf(header2,basesToShow); + coefficients = sysBasesObj.referenceCoefficients; + % To remove numerical noise. + coefficients(abs(coefficients) (up to order 4) for the hexahedron. + + xyzObj = math.RealCoordinates(true); + xyz = xyzObj.coordinates; + x = xyz(1); y = xyz(2); z = xyz(3); + + % We obtain the space as the tensor product of the mixed poly spaces. + sizeBases = 3*order*(order+1)^2; + symPolynomials = sym(zeros(sizeBases,3)); + symPolynomials(1:order*(order+1)^2,1) = math.getMixedPolySpace([x y z], order-1, order, order); + symPolynomials(order*(order+1)^2+1:2*order*(order+1)^2,2) = math.getMixedPolySpace([x y z], order, order-1, order); + symPolynomials(2*order*(order+1)^2+1:end,3) = math.getMixedPolySpace([x y z], order, order, order-1); \ No newline at end of file diff --git a/+bases/showPolynomials.m b/+bases/showPolynomials.m new file mode 100644 index 0000000..4d6bcc5 --- /dev/null +++ b/+bases/showPolynomials.m @@ -0,0 +1,11 @@ +function showPolynomials() + % SHOWBASES displays the order of the polynomials that we follow + % for applying the coefficients. + + symPolynomials = bases.getPolynomials(2); + + for indexPolynomial = 1:size(symPolynomials,1) + fprintf('Polynomial #%2i | [%9s,%9s,%9s]\n',indexPolynomial,symPolynomials(indexPolynomial,:)); + end + +end \ No newline at end of file diff --git a/+geometry/@Edge/Edge.m b/+geometry/@Edge/Edge.m new file mode 100644 index 0000000..6c337cb --- /dev/null +++ b/+geometry/@Edge/Edge.m @@ -0,0 +1,36 @@ +classdef Edge < handle + % This class represents an edge of the mesh defined by two points. They are + % only defined once in the mesh and the elements store references to all + % their edges. + + properties + + % Vertex objects that compose the edge + nodes geometry.Node + + % Local order of the entity. + pOrder double + + % Array of numbers which represent the numbering for FEM matrices (DOF). + dof(:,1) double + + % Unique identifier in the global mesh + id double + + % Identifier of the boundary condition for polyFaces + polyFaceId double + + end + + methods + function edgeObject = Edge(nodes, id) + % Constructor of the edge objects. + + if (nargin > 0) + edgeObject.nodes = nodes; + edgeObject.id = id; + end + end + end + +end \ No newline at end of file diff --git a/+geometry/@Edge/generateDOF.m b/+geometry/@Edge/generateDOF.m new file mode 100644 index 0000000..9928bc9 --- /dev/null +++ b/+geometry/@Edge/generateDOF.m @@ -0,0 +1,7 @@ +function afterDOF = generateDOF(edgeObj, currentDOF) + % GENERATEDOF - returns the last number of dof assigned. + % GENERATEDOF(meshObj) generates two dofs per each edge. + + % For edges, the number of DOF is equal to their order. + edgeObj.dof = currentDOF+1:currentDOF+edgeObj.pOrder; + afterDOF = currentDOF+edgeObj.pOrder; diff --git a/+geometry/@Edge/getJacobian.m b/+geometry/@Edge/getJacobian.m new file mode 100644 index 0000000..dba01c3 --- /dev/null +++ b/+geometry/@Edge/getJacobian.m @@ -0,0 +1,6 @@ +function [jacobianEdge, originCoordinate] = getJacobian(edgeObj) + % GETJACOBIAN returns the mapping between 1D to 3D points for + % the edge passed as parameter. + + jacobianEdge = edgeObj.nodes(2).coordinates - edgeObj.nodes(1).coordinates; + originCoordinate = edgeObj.nodes(1).coordinates; \ No newline at end of file diff --git a/+geometry/@Edge/getLength.m b/+geometry/@Edge/getLength.m new file mode 100644 index 0000000..b9107dd --- /dev/null +++ b/+geometry/@Edge/getLength.m @@ -0,0 +1,11 @@ +function lengthEdge = getLength(edgeObj) + % GETLENGTH gives the physical length for an (array) of edges. + % GETLENGTH(edgeObj) uses the object (or array of objects) as parameter and + % computes the physical length between the coordinates. + + lengthEdge = zeros(1,length(edgeObj)); + for ii = 1:length(edgeObj) + lengthEdge(ii) = norm(edgeObj(ii).nodes(1).coordinates - edgeObj(ii).nodes(2).coordinates); + end + +end diff --git a/+geometry/@Edge/getUnitVector.m b/+geometry/@Edge/getUnitVector.m new file mode 100644 index 0000000..7c62e0a --- /dev/null +++ b/+geometry/@Edge/getUnitVector.m @@ -0,0 +1,13 @@ +function unitVector = getUnitVector(edgeObj) + % GETLENGTH gives the unit vector for an (array) of edges. + % GETLENGTH(edgeObj) uses the object (or array of objects) as parameter and + % computes the unit vector from the coordinates. + + numDimensions = length(edgeObj(1).nodes(1).coordinates); + unitVector = zeros(numDimensions, length(edgeObj)); + for ii = 1:length(edgeObj) + unitVector(:,ii) = edgeObj(ii).nodes(2).coordinates - edgeObj(ii).nodes(1).coordinates; + unitVector(:,ii) = unitVector(:,ii)/norm(unitVector(:,ii)); + end + +end diff --git a/+geometry/@Element/Element.m b/+geometry/@Element/Element.m new file mode 100644 index 0000000..e5517e0 --- /dev/null +++ b/+geometry/@Element/Element.m @@ -0,0 +1,111 @@ +classdef Element < handle + % Class which represents the element in a mesh. It stores the reference to + % all the entities (nodes, edges, faces) which constitute the element. + + properties + + % Nodes which compose the element + nodes geometry.Node + % Edges which compose the element + edges geometry.Edge + % Faces which compose the element + faces geometry.Face + + % Array of interior dofs. + volumeDof double + % Order of the interior element. + pOrder double + % Identifier of the solid (material) which this element belongs to. + solidId = 0; + % Unique identifier in the mesh. + id double; + + end + + methods + function elemObj = Element(varargin) + % Constructor of the elements. Actually, we use it only for + % the reference element since the mesh elements are already + % stored in meshesHalfFilledTAP2021. + + import geometry.*; + + if (nargin > 0) + if (strcmp(varargin{1},'referenceElement')) + % Coordinates of the reference element. + coordinates = [0 0 0; 1 0 0; 1 1 0; 0 1 0; 0 0 1; 1 0 1; 1 1 1; 0 1 1]; + numNodes = 8; + numEdges = 12; + numFaces = 6; + + % Initialization. + objNodes(numNodes) = Node; + objEdges(numEdges) = Edge; + objFaces(numFaces) = Face; + + % Creation of the nodes from the coordinates. + for ii = 1:numNodes + objNodes(ii) = Node(coordinates(ii,:),ii); + end + % We follow the ordering given by getEdgeByLocalNodes. + for ii = 1:numEdges + edgeByLocalNodes = Element.getEdgeByLocalNodes(ii); + objEdges(ii) = Edge(objNodes(edgeByLocalNodes), ii); + end + % We follow the ordering given by getFaceByLocalEdges. + for ii = 1:numFaces + faceByLocalEdges = Element.getFaceByLocalEdges(ii); + objFaces(ii) = Face(objEdges(faceByLocalEdges), ii); + end + % Assignation + elemObj.nodes = objNodes; + elemObj.edges = objEdges; + elemObj.faces = objFaces; + end + end + end + end + + methods (Static) + function edgeByLocalNodes = getEdgeByLocalNodes(edgeNumber) + % GETEDGEBYLOCALNODES This method gives the relation of each edge in + % the element with the nodes of the element. + % GETEDGEBYLOCALNODES(edgeNumber) only returns the nodes + % related to the edge specified in the parameter list. + + edgesByLocalNodes = zeros(12,2); + edgesByLocalNodes(1,:) = [1 2]; + edgesByLocalNodes(2,:) = [2 3]; + edgesByLocalNodes(3,:) = [4 3]; + edgesByLocalNodes(4,:) = [1 4]; + edgesByLocalNodes(5,:) = [5 6]; + edgesByLocalNodes(6,:) = [6 7]; + edgesByLocalNodes(7,:) = [8 7]; + edgesByLocalNodes(8,:) = [5 8]; + edgesByLocalNodes(9,:) = [1 5]; + edgesByLocalNodes(10,:) = [2 6]; + edgesByLocalNodes(11,:) = [3 7]; + edgesByLocalNodes(12,:) = [4 8]; + + edgeByLocalNodes = edgesByLocalNodes(edgeNumber,:); + end + function faceByLocalEdges = getFaceByLocalEdges(faceNumber) + % GETFACEBYLOCALEDGES This method returns the relation of each face + % in the element with the edges of the element. + % GETFACEBYLOCALEDGES(faceNumber) only returns the edges + % related to the face specified in the parameter list. + + facesByLocalEdges = zeros(6,4); + facesByLocalEdges(1,:) = [1 2 3 4]; + facesByLocalEdges(2,:) = [5 6 7 8]; + facesByLocalEdges(3,:) = [1 10 5 9]; + facesByLocalEdges(4,:) = [2 11 6 10]; + facesByLocalEdges(5,:) = [3 11 7 12]; + facesByLocalEdges(6,:) = [4 12 8 9]; + + faceByLocalEdges = facesByLocalEdges(faceNumber,:); + end + end + + +end \ No newline at end of file diff --git a/+geometry/@Element/computeJacobianMatrix.m b/+geometry/@Element/computeJacobianMatrix.m new file mode 100644 index 0000000..f7ec612 --- /dev/null +++ b/+geometry/@Element/computeJacobianMatrix.m @@ -0,0 +1,13 @@ +function jacobMatrix = computeJacobianMatrix(elemObj) + % COMPUTEJACOBIANMATRIX returns the jacobian matrix for straight hexahedra. + % COMPUTEJACOBIANMATRIX(elemObj) uses the object passed as parameter and + % gets the jacobian matrix that represents the mapping between the + % reference element defined in the constructor of element and + % the real element passed as parameter. + + verticesElement = elemObj.getAllCoordinates(); + jacobMatrix = [ verticesElement(2,:) - verticesElement(1,:) + verticesElement(4,:) - verticesElement(1,:) + verticesElement(5,:) - verticesElement(1,:)]; + +end \ No newline at end of file diff --git a/+geometry/@Element/generateVolumeDOF.m b/+geometry/@Element/generateVolumeDOF.m new file mode 100644 index 0000000..dbe836b --- /dev/null +++ b/+geometry/@Element/generateVolumeDOF.m @@ -0,0 +1,8 @@ +function afterDOF = generateVolumeDOF(elemObj, currentDOF) + % GENERATEVOLUMEDOF - returns the last number of dof assigned. + % GENERATEVOLUMEDOF(meshObj) generates six dofs per element. + + sizeDOF = 3*elemObj.pOrder*(elemObj.pOrder-1)^2; + + elemObj.volumeDof = currentDOF+1:currentDOF+sizeDOF; + afterDOF = currentDOF+sizeDOF; diff --git a/+geometry/@Element/getAllCoordinates.m b/+geometry/@Element/getAllCoordinates.m new file mode 100644 index 0000000..4a99be1 --- /dev/null +++ b/+geometry/@Element/getAllCoordinates.m @@ -0,0 +1,11 @@ +function allCoordinates = getAllCoordinates(elemObj) + % GETALLCOORDINATES gets the physical coordinates of all the nodes that + % compose the element. + % GETALLCOORDINATES(elemObj) returns the physical coordinates of all the + % nodes. + + allCoordinates = zeros(8,3); + + for ii = 1:8 + allCoordinates(ii,:) = elemObj.nodes(ii).coordinates; + end \ No newline at end of file diff --git a/+geometry/@Element/getDOF.m b/+geometry/@Element/getDOF.m new file mode 100644 index 0000000..61a9acb --- /dev/null +++ b/+geometry/@Element/getDOF.m @@ -0,0 +1,20 @@ +function dof = getDOF(elemObj) + % GETDOF returns all the dofs within an element + % GETDOF(elemObj) takes the object as parameter and gets the + % list of dofs for edges, faces and volume following the + % order within the element. + + dof = zeros (54,1); + indexDOF = 0; + for indexEdge = 1:length(elemObj.edges) + dof(indexDOF+1:indexDOF+length(elemObj.edges(indexEdge).dof)) = elemObj.edges(indexEdge).dof; + indexDOF = indexDOF + length(elemObj.edges(indexEdge).dof); + end + for indexFace = 1:length(elemObj.faces) + dof(indexDOF+1:indexDOF+length(elemObj.faces(indexFace).dof)) = elemObj.faces(indexFace).dof; + indexDOF = indexDOF + length(elemObj.faces(indexFace).dof); + end + dof (indexDOF+1:indexDOF+length(elemObj.volumeDof)) = elemObj.volumeDof; + indexDOF = indexDOF+length(elemObj.volumeDof); + + dof = dof(1:indexDOF); \ No newline at end of file diff --git a/+geometry/@Face/Face.m b/+geometry/@Face/Face.m new file mode 100644 index 0000000..d26d7be --- /dev/null +++ b/+geometry/@Face/Face.m @@ -0,0 +1,36 @@ +classdef Face < handle + % This class is used to represent a face in the mesh. They are only defined + % once in the mesh and the elements store references to all their faces. + + properties + + % Handles to the edge objects which compose the face. + edges(:,1) geometry.Edge + + % Local order of the entity. + pOrder double + + % Array of numbers which represent the numbering for + % FEM matrices (DOF). + dof(:,1) double + + % Identifier in the global mesh + id double + + % Identifier of the boundary condition for polyFaces + polyFaceId double + + end + + methods + function faceObject = Face(edges,id) + % Constructor of the face objects. + + if (nargin > 0) + faceObject.edges = edges; + faceObject.id = id; + end + end + end + +end \ No newline at end of file diff --git a/+geometry/@Face/generateDOF.m b/+geometry/@Face/generateDOF.m new file mode 100644 index 0000000..4a68ed7 --- /dev/null +++ b/+geometry/@Face/generateDOF.m @@ -0,0 +1,7 @@ +function afterDOF = generateDOF(faceObj, currentDOF) + % GENERATEDOF - returns the last number of dof assigned. + % GENERATEDOF(meshObj) generates four dofs per each face. + + sizeDOF = 2*faceObj.pOrder*(faceObj.pOrder-1); + faceObj.dof = currentDOF+1:currentDOF+sizeDOF; + afterDOF = currentDOF+sizeDOF; diff --git a/+geometry/@Mesh/Mesh.m b/+geometry/@Mesh/Mesh.m new file mode 100644 index 0000000..760cdd7 --- /dev/null +++ b/+geometry/@Mesh/Mesh.m @@ -0,0 +1,35 @@ +classdef Mesh < handle + % This is a class used to store all the references to nodes, edges, faces, elements, + % materials, and boundary conditions presents in a mesh. + + properties + + % Reference to nodes. + nodes(:,1) geometry.Node + % Reference for edges. + edges(:,1) geometry.Edge + % Reference for faces. + faces(:,1) geometry.Face + % Reference for elements. + elements(:,1) geometry.Element + % Reference for the solids in the mesh. + solids(:,1) geometry.Solid + % Reference for the boundary conditions in the mesh. + polyFaces(:,1) geometry.PolyFace + + % Number of DOF present in the mesh. + numDOF double + % Settings of the problem where mesh is the domain. + settingsObj settings.SystemSettings + + end + + methods + function meshObj = Mesh() + % This constructor is the default one since the mesh are loaded from + % memory. + + end + end + +end \ No newline at end of file diff --git a/+geometry/@Mesh/generateDOF.m b/+geometry/@Mesh/generateDOF.m new file mode 100644 index 0000000..d74c735 --- /dev/null +++ b/+geometry/@Mesh/generateDOF.m @@ -0,0 +1,21 @@ +function currentDOF = generateDOF(meshObj) + % GENERATEDOF - returns the final number of DOF that are stored + % inside the mesh. + % GENERATEDOF(meshObj) runs through all the entities (edges, faces and elements) + % in the mesh to provide a unique identifier for each dof in each entity. + + currentDOF = 0; + + for ii = 1:length(meshObj.edges) + currentDOF = meshObj.edges(ii).generateDOF(currentDOF); + end + + for ii = 1:length(meshObj.faces) + currentDOF = meshObj.faces(ii).generateDOF(currentDOF); + end + + for ii = 1:length(meshObj.elements) + currentDOF = meshObj.elements(ii).generateVolumeDOF(currentDOF); + end + + meshObj.numDOF = currentDOF; diff --git a/+geometry/@Mesh/getLengthStatistics.m b/+geometry/@Mesh/getLengthStatistics.m new file mode 100644 index 0000000..56eb5ff --- /dev/null +++ b/+geometry/@Mesh/getLengthStatistics.m @@ -0,0 +1,11 @@ +function [avgLength, maxLength, minLength] = getLengthStatistics(meshObj) + % GETLENGTHSTATISTICS gets the average, maximum and minimum lengths of the + % edges in the mesh. + % GETLENGTHSTATISTICS(meshObj) runs through all the edges in the mesh passed + % as parameter and returns the average, maximum and minimum lengths. + + lengthEdges = meshObj.edges.getLength(); + + avgLength = mean(lengthEdges); + maxLength = max(lengthEdges); + minLength = min(lengthEdges); \ No newline at end of file diff --git a/+geometry/@Mesh/getNNZ.m b/+geometry/@Mesh/getNNZ.m new file mode 100644 index 0000000..5843449 --- /dev/null +++ b/+geometry/@Mesh/getNNZ.m @@ -0,0 +1,15 @@ +function nonZeros = getNNZ(meshObj) + % GETNNZ returns an upper bound for the number of non-zeros + % of a given problem. + % GETNNZ(meshObj) takes the mesh and runs through all the elements + % counting the number of dofs associated to the element. + + nonZerosArray = zeros(meshObj.numDOF,1); + + for ii = 1:length(meshObj.elements) + dofElement = meshObj.elements(ii).getDOF(); + nonZerosArray (dofElement) = nonZerosArray(dofElement) + length(dofElement); + end + nonZeros = sum(nonZerosArray); + + diff --git a/+geometry/@Node/Node.m b/+geometry/@Node/Node.m new file mode 100644 index 0000000..d1de38f --- /dev/null +++ b/+geometry/@Node/Node.m @@ -0,0 +1,27 @@ +classdef Node < handle + % Class used to store the coordinates and the identifiers of the node in the mesh. + + properties + % Physical coordinates for the point. + % It is a double array of 3 dimensions. + coordinates(:,1) double + + % Identifier in the global mesh + id double + + % Boundary condition from polyFaces + polyFaceId double + end + + methods + function nodeObject = Node(coordinates, id) + % Constructor of the node objects. + if (nargin > 0) + nodeObject.coordinates = coordinates; + nodeObject.id = id; + nodeObject.polyFaceId = 0; + end + end + end + +end \ No newline at end of file diff --git a/+geometry/@PolyFace/PolyFace.m b/+geometry/@PolyFace/PolyFace.m new file mode 100644 index 0000000..f0ec707 --- /dev/null +++ b/+geometry/@PolyFace/PolyFace.m @@ -0,0 +1,24 @@ +classdef PolyFace < handle + % Class used to store the boundary conditions present in the problem. + + properties + % Identifier for the PolyFace + id double + % Name for the PolyFace. + name + % String which define the type of the PolyFace + polyFaceType + % Identifiers for the faces which belong to the PolyFace. + faces(:,1) double + + end + + methods + function polyFaceObject = PolyFace() + % This constructor is the default one since the mesh (and poly faces + % associated) are loaded from memory. + + end + end + +end \ No newline at end of file diff --git a/+geometry/@Solid/Solid.m b/+geometry/@Solid/Solid.m new file mode 100644 index 0000000..7811339 --- /dev/null +++ b/+geometry/@Solid/Solid.m @@ -0,0 +1,20 @@ +classdef Solid < handle + % Class used to store the solids (materials) present in the problem. + + properties + % Name for the Solid. + name + % Identifiers for the elements (or faces in 2D) which belong to the Solid. + elements(:,1) double + % Material of the solid + materialObj settings.Material + end + + methods + function solidObject = Solid() + % This constructor is the default one since the mesh (and solids + % associated) are loaded from memory. + end + end + +end \ No newline at end of file diff --git a/+math/@Polynomial/Polynomial.m b/+math/@Polynomial/Polynomial.m new file mode 100644 index 0000000..d3b0902 --- /dev/null +++ b/+math/@Polynomial/Polynomial.m @@ -0,0 +1,59 @@ +classdef Polynomial < handle + + % Class used to represent polynomials as a multivariate array + % of coefficients to speed up the computations. + % This coefficients may be symbolic (very slow) or double numbers. + + properties + % Double array with different size depending on the variables involved + % and the maximum order of the polynomial. + coefficients + end + + methods + function polynomialObject = Polynomial(maxOrder, realExpression) + % maxOrder is an array which specifies the maximum order for all the + % polynomials involved. + % [coordinates] and [symExpression] are optional, and when included + % the polynomial is created with this coefficients using poly2list. + % deBoorFlag is a flag which should be used when the polynomials are stored. + if (nargin > 0) + % Array of three dimensions (typically) + if (length(maxOrder) == 1) + polynomialObject.coefficients = zeros(maxOrder+1,1); + % polynomialObject.symCoefficients = sym(zeros(maxOrder+1,1)); + else + polynomialObject.coefficients = zeros(maxOrder+1); + % polynomialObject.symCoefficients = sym(zeros(maxOrder+1)); + end + + if ((length(maxOrder) < 1) || (length(maxOrder) > 3)) + error ('The order of the polynomial representation have to be specified for each dimension'); + end + + coordsObj = math.RealCoordinates(true); + for ii = 1:length(coordsObj.coordinates) + assume(coordsObj.coordinates(ii), 'real'); + end + + if (~isequaln(simplify(realExpression),sym(0))) + coords = coordsObj.coordinates; + + [coefficients, exponents] = coeffs(realExpression); + % Get place of coefficients. + for ii = 1:length(coefficients) + % Working from R2018a + polynomialObject.coefficients( ... + polynomialDegree(exponents(ii),coords(1))+1,... + polynomialDegree(exponents(ii),coords(2))+1,... + polynomialDegree(exponents(ii),coords(3))+1) = ... + double(coefficients(ii)); + end + + end + + end + end + end + +end \ No newline at end of file diff --git a/+math/@Polynomial/evaluatePolynomial.m b/+math/@Polynomial/evaluatePolynomial.m new file mode 100644 index 0000000..427b235 --- /dev/null +++ b/+math/@Polynomial/evaluatePolynomial.m @@ -0,0 +1,27 @@ +function evalPoly = evaluatePolynomial(polyObj, r) + % EVALUATEPOLYNOMIAL gives the evaluation of the polynomial + % object in a point r. + % EVALUATEPOLYNOMIAL(polyObj, r) runs through all the coefficients in the + % polynomial and use the vector coordinate r to return the value of the polynomial. + + % To speed up computations. + tolerance = 1e-16; + + evalPoly = 0; + coefficients = polyObj.coefficients; + + lengthT1 = size(coefficients,1); + lengthT2 = size(coefficients,2); + lengthT3 = size(coefficients,3); + for t1Index = 1:lengthT1 + for t2Index = 1:lengthT2 + for t3Index = 1:lengthT3 + if (abs(coefficients(t1Index,t2Index,t3Index))>tolerance) + evalPoly = evalPoly + coefficients(t1Index,t2Index,t3Index)*... + r(1)^(t1Index-1)*r(2)^(t2Index-1)*r(3)^(t3Index-1); + end + end + end + end + + end \ No newline at end of file diff --git a/+math/@Polynomial/getSymPolynomial.m b/+math/@Polynomial/getSymPolynomial.m new file mode 100644 index 0000000..9e47051 --- /dev/null +++ b/+math/@Polynomial/getSymPolynomial.m @@ -0,0 +1,29 @@ +function symPoly = getSymPolynomial(polyObj) + % GETSYMPOLYNOMIAL returns a symbolic representation of the polynomial + % object. + % GETSYMPOLYNOMIAL(polyObj) use symbolic real coordinates to return the + % symbolic representation of the polynomial (useful for debugging purposes). + % The polynomial object may be an array. + + % Initialization. + symPoly = sym(zeros(size(polyObj))); + + % Run through all the polynomial objects passed as parameter. + for indexPolyRow = 1:size(polyObj,1) + for indexPolyCol = 1:size(polyObj,2) + realObj = math.RealCoordinates(); + realCoords = realObj.coordinates; + + coefficients = polyObj(indexPolyRow,indexPolyCol).coefficients; + + % Run through all the coefficients. + for t1Index = 1:size(coefficients,1) + for t2Index = 1:size(coefficients,2) + for t3Index = 1:size(coefficients,3) + symPoly(indexPolyRow, indexPolyCol) = symPoly(indexPolyRow, indexPolyCol) + coefficients(t1Index,t2Index,t3Index)*... + realCoords(1)^(t1Index-1)*realCoords(2)^(t2Index-1)*realCoords(3)^(t3Index-1); + end + end + end + end + end \ No newline at end of file diff --git a/+math/@RealCoordinates/RealCoordinates.m b/+math/@RealCoordinates/RealCoordinates.m new file mode 100644 index 0000000..d678bd1 --- /dev/null +++ b/+math/@RealCoordinates/RealCoordinates.m @@ -0,0 +1,58 @@ +classdef RealCoordinates < handle + % This class is used to store the symbolic representation of x, y, and z. + + properties + % Symbolic array with different size depending on the geometry of the object. + coordinates sym + end + + methods + function realCoordinatesObject = RealCoordinates(loadSymVariables) + % Constructor which defines the array of symbolic variables x,y,z + % used in the code, and if they exist, they load them (it is a little + % slow for MATLAB). + + if (nargin > 0) + loadedSym = false; + if (loadSymVariables) + realName = 'realCoordinates'; + fullFile = strcat('+math/@RealCoordinates/precomputedReal/',realName,'.mat'); + if (isfile(fullFile)) + load(fullFile,'tempReal'); + realCoordinatesObject.coordinates = tempReal.coordinates; + loadedSym = true; + end + end + + if (~loadedSym) + + saveFunctions = true; + + syms x y z + assume(x, 'real'); + assume(y, 'real'); + assume(z, 'real'); + realCoordinatesObject.coordinates = [x y z]; + + if (saveFunctions) + realCoordinatesObject.saveRealCoordinates(); + end + end + end + end + + function saveRealCoordinates (realCoordinatesObject) + % Function used to store the symbolic variables used in RealCoordinates. + + realName = 'realCoordinates'; + tempReal = 'tempReal'; + functionsStruct.(tempReal) = realCoordinatesObject; + fullFile = fullfile('+math/@RealCoordinates/precomputedReal/',realName); + if (~exist('+math/@RealCoordinates/precomputedReal/', 'dir')) + mkdir ('+math/@RealCoordinates/precomputedReal/'); + end + save (fullFile, '-struct', 'functionsStruct'); + end + end + +end \ No newline at end of file diff --git a/+math/getFullPolySpace.m b/+math/getFullPolySpace.m new file mode 100644 index 0000000..4f499d0 --- /dev/null +++ b/+math/getFullPolySpace.m @@ -0,0 +1,11 @@ +function polySpace = getFullPolySpace(t, order) + % GETFULLPOLYSPACE Returns Q(order) for all the coordinates. + % GETFULLPOLYSPACE(t, order) returns the polynomial space for each component + % up to 3 dimensions. + % t is a 1D symbolic variable. + % See also GETMIXEDPOLYSPACE + + polySpace = sym(zeros(1,(order+1))); + for pp = 0:order + polySpace(pp+1) = t^pp; + end diff --git a/+math/getGaussPoints.m b/+math/getGaussPoints.m new file mode 100644 index 0000000..a5fc196 --- /dev/null +++ b/+math/getGaussPoints.m @@ -0,0 +1,29 @@ +function [weights,points] = getGaussPoints(order) + % GETGAUSSPOINT gives the weights and 3D points to perform + % numerical integration on one hexahedron. + % GETGAUSSPOINT(order) returns the points that we need to do + % a numerical integration for a function up to the order passed as + % parameter. These points are obtained as the tensor product of the + % 1D gauss points for the segment and the 2D points for the quadrilateral. + % See also GETGAUSSPOINTFORSEGMENT, GETGAUSSPOINTSFORQUADRILATERAL + + + [weightsSegment, pointsSegment] = math.getGaussPointsForSegment(order); + [weights2D, points2D] = math.getGaussPointsForQuadrilateral(order); + + num1DGauss = length(weightsSegment); + num2DGauss = length(weights2D); + + weights = zeros(num1DGauss*num2DGauss,1); + points = zeros(num1DGauss*num2DGauss,3); + + for index_points = 1:num1DGauss + weights((index_points-1)*num2DGauss+1:index_points*num2DGauss) = weightsSegment(index_points)*weights2D(1:num2DGauss); + points((index_points-1)*num2DGauss+1:index_points*num2DGauss,1:2) = points2D(1:num2DGauss,:); + points((index_points-1)*num2DGauss+1:index_points*num2DGauss,3) = pointsSegment(index_points); + end + +end + + + diff --git a/+math/getGaussPointsForQuadrilateral.m b/+math/getGaussPointsForQuadrilateral.m new file mode 100644 index 0000000..11ddf62 --- /dev/null +++ b/+math/getGaussPointsForQuadrilateral.m @@ -0,0 +1,24 @@ +function [weights, points] = getGaussPointsForQuadrilateral(order) + % GETGAUSSPOINTFORQUADRILATERAL gives the weights and 2D points to perform + % numerical integration on one quadrilateral face. + % GETGAUSSPOINTFORQUADRILATERAL(order) returns the points that we need to do + % a numerical integration for a function up to the order passed as + % parameter. These points are obtained as the tensor product of the + % 1D gauss points for the segment. + % See also GETGAUSSPOINTFORSEGMENT + + % We get the corresponding points for one segment. + [weightsSegment, pointsSegment] = math.getGaussPointsForSegment(order); + + num1DGauss = length(weightsSegment); + + weights = zeros(num1DGauss^2,1); + points = zeros(num1DGauss^2,2); + + % We implement here the tensor product of the 1D points. + for indexPoints = 1:num1DGauss + weights((indexPoints-1)*num1DGauss+1:indexPoints*num1DGauss) = weightsSegment(indexPoints)*weightsSegment(1:num1DGauss); + points((indexPoints-1)*num1DGauss+1:indexPoints*num1DGauss,1) = pointsSegment(1:num1DGauss); + points((indexPoints-1)*num1DGauss+1:indexPoints*num1DGauss,2) = pointsSegment(indexPoints); + end +end \ No newline at end of file diff --git a/+math/getGaussPointsForSegment.m b/+math/getGaussPointsForSegment.m new file mode 100644 index 0000000..0ffedb5 --- /dev/null +++ b/+math/getGaussPointsForSegment.m @@ -0,0 +1,32 @@ +function [weights, points] = getGaussPointsForSegment(order) + % GETGAUSSPOINTFORSEGMENT gives the weights and 1D points to perform + % numerical integration at one segment. + % GETGAUSSPOINTFORSEGMENT(order) returns the points that + % we need to do a numerical integration for a function + % up to the order passed as parameter. + + % In terms of the order we return a different number of points. + switch order + + case num2cell([0 1]) + weights = 1; + points = 0.555555555555555; + + case num2cell([2 3]) + weights(1:2) = 1/2; + points(1) = 0.788675134594813; + points(2) = 0.211324865405187; + + case num2cell([4 5]) + weights = [ 0.277777777777778 + 0.444444444444444 + 0.277777777777778]; + + points = [ 0.112701665379259 + 0.500000000000000 + 0.887298334620741]; + otherwise + error("Order not coded") + end + +end \ No newline at end of file diff --git a/+math/getMixedPolySpace.m b/+math/getMixedPolySpace.m new file mode 100644 index 0000000..7c3fe6b --- /dev/null +++ b/+math/getMixedPolySpace.m @@ -0,0 +1,46 @@ +function polySpace = getMixedPolySpace(t,order1,order2,order3) + % GETMIXEDPOLYSPACE Returns Q(order1,order2,order3) for t. + % GETMIXEDPOLYSPACE(t,order1,order2,order3) returns the mixed polynomial + % space of Q(order1,order2,order3). + % t is a 3D array of symbolic variables + % order1 is the order of the first coordinate + % order2 is the order of the second coordinate + % order3 is the order of the third coordinate + % + % See also GETFULLPOLYSPACE + + if ((min([order1, order2, order3])<0) && (max([order1, order2, order3])>3)) + disp('This order is not supported'); + return + end + + % We get the full poly space for all the components. + polySpace = sym(zeros((order1+1)*(order2+1)*(order3+1),1)); + polySpace1 = math.getFullPolySpace(t(1),order1); + polySpace2 = math.getFullPolySpace(t(2),order2); + polySpace3 = math.getFullPolySpace(t(3),order3); + + + % We obtain the full tensor product of all the spaces + mixedProducts_t1t2 = polySpace1.'*polySpace2; + mixedProducts = sym(zeros(size(mixedProducts_t1t2,1),size(mixedProducts_t1t2,2),length(polySpace3))); + for ii = 1:length(polySpace3) + mixedProducts(:,:,ii) = polySpace3(ii)*mixedProducts_t1t2; + end + + % The order of the resulting matrix is the index. So, only products with + % order less than the actual order requested (variable called here "order") + % will be stored. Note that -1 has to be substracted to all the indices + % since the order for the element (1,1) is 0. + indexPolySpace = 1; + for ii = 1:length(polySpace3) + for jj = 1:length(polySpace2) + for kk = 1:length(polySpace1) + % We remove the higher-order terms. + if ((ii+jj+kk-3)<=(order1+order2+order3)) + polySpace(indexPolySpace) = mixedProducts(kk,jj,ii); + indexPolySpace = indexPolySpace + 1; + end + end + end + end \ No newline at end of file diff --git a/+math/transform2Dto3DInReferenceElement.m b/+math/transform2Dto3DInReferenceElement.m new file mode 100644 index 0000000..664b780 --- /dev/null +++ b/+math/transform2Dto3DInReferenceElement.m @@ -0,0 +1,25 @@ +function points3D = transform2Dto3DInReferenceElement(points2D, faceNumber) + % TRANSFORM2DTO3DINREFERENCEELEMENT - Transforms 2D points defined on some + % face of the hexahedron into its 3D representation. + % + % TRANSFORM2DTO3DINREFERENCEELEMENT(points2D, faceNumber) takes the 2D representation + % of the point and returns the corresponding 3D representation for a given face + % specified by faceNumber. The same numbering as in GETFACEBYLOCALEDGES is used. + + + switch faceNumber + case 1 + points3D = [points2D(1) points2D(2) 0]; + case 2 + points3D = [points2D(1) points2D(2) 1]; + case 3 + points3D = [points2D(1) 0 points2D(2)]; + case 4 + points3D = [1 points2D(1) points2D(2)]; + case 5 + points3D = [points2D(1) 1 points2D(2)]; + case 6 + points3D = [0 points2D(1) points2D(2)]; + otherwise + error('A hexahedron only has 6 faces.') + end diff --git a/+settings/@Material/Material.m b/+settings/@Material/Material.m new file mode 100644 index 0000000..ed822a2 --- /dev/null +++ b/+settings/@Material/Material.m @@ -0,0 +1,26 @@ +classdef Material < handle + + % Class to represent the material of the electromagnetic problems. + % Inhomogeneous and anisotropic material is supported (is stored in + % the parameters and called with evaluateMaterial). + + properties + % Name for the material. + name char + % Reference to the material property object. + epsilonRelative settings.MaterialProperty + % Reference to the material property object. + muRelative settings.MaterialProperty + % Reference to the material property object. + sigma settings.MaterialProperty + + end + + methods + function materialObject = Material() + % Default constructor since everything is loaded from the .mat + % container. + end + end + +end \ No newline at end of file diff --git a/+settings/@MaterialProperty/MaterialProperty.m b/+settings/@MaterialProperty/MaterialProperty.m new file mode 100644 index 0000000..ec89f6a --- /dev/null +++ b/+settings/@MaterialProperty/MaterialProperty.m @@ -0,0 +1,19 @@ +classdef MaterialProperty < handle + + % Class to represent the material of the electromagnetic problems. + % Inhomogeneous and anisotropic material is supported (is stored in + % the parameters and called with evaluateProperty). + + properties + % This might be a double, a symbolic expression or a math.Polynomial. + property + end + + methods + function materialPropertyObject = MaterialProperty() + % Default constructor since everything is loaded from the .mat + % container. + end + end + +end \ No newline at end of file diff --git a/+settings/@MaterialProperty/evaluateProperty.m b/+settings/@MaterialProperty/evaluateProperty.m new file mode 100644 index 0000000..a779423 --- /dev/null +++ b/+settings/@MaterialProperty/evaluateProperty.m @@ -0,0 +1,4 @@ +function propertyEvaluated = evaluateProperty (materialPropertyObject) + % EVALUATEPROPERTY returns the value of the material property. + + propertyEvaluated = materialPropertyObject.property; \ No newline at end of file diff --git a/+settings/@SystemSettings/SystemSettings.m b/+settings/@SystemSettings/SystemSettings.m new file mode 100644 index 0000000..5dd4dd7 --- /dev/null +++ b/+settings/@SystemSettings/SystemSettings.m @@ -0,0 +1,21 @@ +classdef SystemSettings < handle + + % Class to represent settings of the problem (order of the basis functions, materials, frequency and curved tailored surfaces if present) + + properties + % Order present in the mesh + % This would be an array of #elements x #order for each entity if necessary. + % Matlab is not really good at dealing with different types of variables so let's restrict to double + pOrder double + % Electromagnetic parameters of the materials + materials settings.Material + end + + methods + function systemSettingsObject = SystemSettings() + % Default constructor since everything is loaded from the .mat + % container. + end + end + +end \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..c93287f --- /dev/null +++ b/README.md @@ -0,0 +1,8 @@ +Demonstration of Second-Order Nédélec Curl-Conforming Hexahedral element +=== + +The main executable is `mainHalfFilled.m`, that generates from scratch the +basis functions and use them to solve the generalized eigenvalue problem +used for verification in the numerical results (specifically, a half-filled +cavity). This script also shows the coefficients for the different basis +shown in Table II. \ No newline at end of file diff --git a/getKcFromHalfFilledCavity.m b/getKcFromHalfFilledCavity.m new file mode 100644 index 0000000..f173ba9 --- /dev/null +++ b/getKcFromHalfFilledCavity.m @@ -0,0 +1,39 @@ +function kc = getKcFromHalfFilledCavity(a, b, c, h, epr, sigmaCond, m, n, axisLimits, resolution) + % GETKCFROMHALFFILLEDCAVITY Taken from https://doi.org/10.5281/zenodo.3885581 + + kc = 0; + + syms kcsquared kcsym + + % Lossy material + if (sigmaCond > eps) + lhsEquation = sqrt(kcsym^2-(m*pi/a)^2-(n*pi/b)^2)*cot(sqrt(kcsym^2-(m*pi/a)^2-(n*pi/b)^2)*(c-h)); + eprMod = epr - 1i*sigmaCond*4*pi*1e-7*299792458/kcsym; + rhsEquation = -sqrt(eprMod*kcsym^2-(m*pi/a)^2-(n*pi/b)^2)*cot(sqrt(eprMod*kcsym^2-(m*pi/a)^2-(n*pi/b)^2)*h); + else + lhsEquation = sqrt(kcsquared-(m*pi/a)^2-(n*pi/b)^2)*cot(sqrt(kcsquared-(m*pi/a)^2-(n*pi/b)^2)*(c-h)); + eprMod = epr; + rhsEquation = -sqrt(eprMod*kcsquared-(m*pi/a)^2-(n*pi/b)^2)*cot(sqrt(eprMod*kcsquared-(m*pi/a)^2-(n*pi/b)^2)*h); + end + + if (nargin > 9) + + digitsOld = digits(resolution); + equation = lhsEquation-rhsEquation == 0; + if (sigmaCond > eps) + kc = (vpasolve(equation, kcsym, axisLimits)); + else + kc = (sqrt(vpasolve(equation, kcsquared, axisLimits.^2))); + end + digits(digitsOld); + else + close all + fplot(real(lhsEquation),axisLimits) + hold on + fplot(real(rhsEquation),axisLimits); + fplot(imag(rhsEquation),axisLimits); + xlabel('k_c') + ylabel('Equation') + legend('LHS','RHS','Location','best') + set(gca, 'fontsize', 22) + end \ No newline at end of file diff --git a/mainHalfFilled.m b/mainHalfFilled.m new file mode 100644 index 0000000..abd1dff --- /dev/null +++ b/mainHalfFilled.m @@ -0,0 +1,78 @@ +%% Convergence for a half-filled cavity. +% This script is used as a demonstration of the capabilities of the basis functions +% in the submission. +% Results from fig. 5 for the half-filled cavity are obtained running this script. + +%% Definition of parameters. +% Parameters for the half-filled cavity as defined in https://doi.org/10.1002/mmce.22342 +a = 0.01; +b = 0.001; +c = 0.01; +h = 0.005; +epr = 2; +sigmaCond = 0; + +%% Obtention of analytic eigenvalues. +% Eigenvalues for the first six modes. +k101 = getKcFromHalfFilledCavity(a, b, c, h, epr, sigmaCond, 1, 0, 300, 16); +k102 = getKcFromHalfFilledCavity(a, b, c, h, epr, sigmaCond, 1, 0, 600, 16); +k103 = getKcFromHalfFilledCavity(a, b, c, h, epr, sigmaCond, 1, 0, 800, 16); +k201 = getKcFromHalfFilledCavity(a, b, c, h, epr, sigmaCond, 2, 0, 540, 16); +k202 = getKcFromHalfFilledCavity(a, b, c, h, epr, sigmaCond, 2, 0, 760, 16); +k301 = getKcFromHalfFilledCavity(a, b, c, h, epr, sigmaCond, 3, 0, 750, 16); +numberKc = 6; +kc = [k101 k201 k102 k301 k202 k103]; + +%% Loading of the hexahedral meshes. +% In this container are included the meshes used for the results. +load('meshesHalfFilledTAP2021.mat','projectMeshes'); +lengthMeshes = length(projectMeshes); +% Container to store the error for the different eigenvalues. +errorHalfConvergence = zeros(numberKc,lengthMeshes); +% Container for the average size of the elements (x-axis in fig. 5). +lengthHalfConvergence = zeros(3,lengthMeshes); + +%% Obtention of the basis functions. +% For demonstration purposes. +disp("Removing the bases stored...") +bases.SystematicBases.removeBases(); +% The bases are defined as an object, and here we define the polynomial q in the +% dofs (6), (7), and (8). +sysBasesObj = bases.SystematicBases(); +disp("Generating the new coefficients...") +% We obtain the coefficients as the dual basis to the Nédélec dofs. +sysBasesObj.getCoefficients(); +% We store the coefficients for using them later. +disp("Storing the coefficients generated...") +sysBasesObj.saveCoefficients(); +disp("Showing the coefficients of the polynomials...") +bases.showPolynomials(); +sysBasesObj.showBases(); + +%% Running through all the meshes. +for indexMesh = 1:length(projectMeshes) + + meshObj = projectMeshes(indexMesh); + + % We generate the dofs inside the mesh. + numDOF = meshObj.generateDOF(); + disp(['Numbering of DOF finished: ', num2str(numDOF), ' unknowns.']) + + % We obtain the global stiffness and mass matrices. + FEMatrices = assembly.getGlobalMatrices(meshObj); + globalStiff = FEMatrices.stiffness; globalMass = FEMatrices.mass; + + % We provide an initial guesh for the sparse eigensolver. + sigmaEstimation = double(kc(1)^2+1.5*(kc(end)-kc(1))^2); + eigsNum = sort(eigs(globalStiff, globalMass, numberKc, sigmaEstimation)); + + % We compute the error for the generalized eigenvalue problem. + errorHalfConvergence (:, indexMesh) = abs((kc.^2)'-eigsNum)./abs((kc.^2)'); + disp(['Mean error obtained: ', num2str(mean(errorHalfConvergence (:, indexMesh)),16)]) + + % We get here some statistics from the mesh to plot the results. + [lengthHalfConvergence(1,indexMesh), lengthHalfConvergence(2,indexMesh), lengthHalfConvergence(3,indexMesh)] = meshObj.getLengthStatistics(); +end + +%% Plotting the results. +plotHalfFilledCavity(errorHalfConvergence,lengthHalfConvergence); diff --git a/meshesHalfFilledTAP2021.mat b/meshesHalfFilledTAP2021.mat new file mode 100644 index 0000000..17f3ade Binary files /dev/null and b/meshesHalfFilledTAP2021.mat differ diff --git a/plotHalfFilledCavity.m b/plotHalfFilledCavity.m new file mode 100644 index 0000000..dff9638 --- /dev/null +++ b/plotHalfFilledCavity.m @@ -0,0 +1,39 @@ +function plotHalfFilledCavity (errorHalfToPlot, lengthHalfToPlot) + % PLOTHALFFILLEDCAVITY generates part of the Fig. 5 + % PLOTHALFFILLEDCAVITY (errorHaltToPlot, lengthHalfToPlot) takes + % as parameter the results from mainHalfFilledCavity and represent + % them in a log-log scale. + % See also MAINHALFFILLED. + + close all + + set(0,'defaultAxesFontName', 'CMU Serif') + + figHandle = figure; + set(figHandle, 'defaultLegendAutoUpdate', 'off') + + loglog(1./lengthHalfToPlot(1,:), errorHalfToPlot(1,:), '-.o', 'linewidth', 2) + hold on + + title('First eigenvalue error for half-filled cavity') + grid on + xlabel('1/h') + ylabel('Relative error') + set(gca,'fontsize',16) + + xlim([1 3000]) + ylim([8e-7 1e-1]) + lengthOriginTriangle = 20; + error1 = 1e-4; + error2 = 1e-6; + + k = 4; + colorOrder4 = [0.4940, 0.1840, 0.5560]; + line([lengthOriginTriangle lengthOriginTriangle], [error1 error2], 'linewidth', 2, 'Color', colorOrder4) + % Pendiente. + length2 = 1/10^(log10(1/lengthOriginTriangle)-(log10(error1)-log10(error2))/k); + h1 = loglog([lengthOriginTriangle length2], [error1 error2]); + set(h1, 'Color', colorOrder4) + set(h1, 'linewidth', 2) + line([lengthOriginTriangle length2], [error2 error2], 'linewidth', 2, 'Color', colorOrder4) + text(lengthOriginTriangle+10, error1/2, ' k=4', 'fontsize', 18, 'fontweight', 'bold', 'fontname','CMU Serif', 'Color', colorOrder4) \ No newline at end of file