Implementación eficiente en MATLAB del Método del Residuo Equilibrado

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Proyecto Fin de Grado Grado en Ingeniería Aeroespacial Implementación eficiente en MATLAB del Método del Residuo Equilibrado. Autor: Antonio Jesús Camúñez Delgado Tutor: Juan Bosco García Archilla Dpto. Matemática Aplicada II Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Proyecto Fin de Grado Grado en Ingeniería Aeroespacial Implementación eficiente en MATLAB del Método del Residuo Equilibrado. Autor: Antonio Jesús Camúñez Delgado Tutor: Juan Bosco García Archilla Catedrático Dpto. Matemática Aplicada II Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Proyecto Fin de Grado: Implementación eficiente en MATLAB del Método del Residuo Equilibrado. Autor: Antonio Jesús Camúñez Delgado Tutor: Juan Bosco García Archilla El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos En primer lugar, a mi tutor por su gran profesionalidad y sinceridad desde el primer días hasta el último. A mi familia por confiar ciegamente en cada una de mis decisiones y ayudarme en todo lo posible. A Susana porque sin su apoyo y su forma de hacerme pensar no podría haber acabado este proyecto Antonio J. Camúñez Delgado Sevilla, 2021 I Resumen El Método de lo Elementos Finitos es una herramienta muy extendida en el campo de la ingeniería para encontrar soluciones aproximadas a problemas con ecuaciones en derivadas parciales, que no suelen tener solución analítica. Hay diversas formar de estimar el error de los resultados calculados, uno de ellos es el Método del Residuo Equilibrado. Tras exponer las bases teóricas, se utilizará MATLAB para elaborar un código que integre el Método del Residuo Equilibrado. Siendo la finalidad del proyecto llegar a elaborar un código definitivo cuyos tiempos tiempos de cálculo sean prácticos, demostrando que se puede conseguir una implementación en MATLAB. Para conseguirlo se estudiarán qué factores ralentizan los cálculos y se irán reduciendo su influencia al mínimo posible. III Introducción En la ingeniería existen muchos métodos analíticos que permiten obtener e forma directa el resultado exacto de un problema. El inconveniente de éstos es que no todos los problemas se pueden resolver se esta forma, ya que, esas soluciones analíticas están asociadas a un tipo de problema con unas dimensiones y unas condiciones de contorno en concreto, como por ejemplo la flecha de una viga empotrada en un extremo. Figura 0.1 Viga en voladizo. En los casos prácticos con los que se encuentra un ingeniero casi siempre hay envueltas ecuaciones en derivadas parciales, con cierta complejidad, aplicadas sobre un amplio abanico de situaciones. Luego, no hay una solución analítica para cada una de ellas. Para obtener una aproximación de la solución se recurre a métodos numéricos, que utilizan el cálculo computacional para obtenerla. De la mano de la solución aproximada nace el concepto de error, que será la diferencia entre la solución numérica y la real. Uno de estos procedimientos más usados sobre muchos campos de la ingeriería es el Método de los Elementos Finitos. Éste se caaracteriza por dividir el cuerpo estudiado en partes más pequeñas, llamadas elementos, dentro de las que se expresa la solución únicamente como función de unos valores en sus vértices que se utilizan para definir funciones, que generalmente serán lineales o cuadráticas. Es decir, la solución aproximada será la combinación de la de cada elemento, teniendo una forma sencilla dentro de cada uno.La válidez de esta aproximación vendrá determinada por el valor del error. Ahí es donde entran el juego los métodos de estimación del error a posteriori, que reutlizan datos calculados durante la obtención de la solución aproximada con elementos finitos. El Método del Residuo Equilibrado es uno de ellos, con el que se obtiene una estimación de la aportación al error de cada elemento gracias a plantear al problema de una forma en concreto. Al conocer la validez de la solución en cada elemento se puede decidir si es necesario dividir algunos en partes más pequeñas, con el fin de llegar a una solución que se 1 acerque lo suficiente a la exacta. El tiempo de cálculo es un factor crucial al usar estas herramientas, estando altamente relacionados con la productividad de un ingeniero en su jornada laboral. De modo que tanto los procesos para obtener la aproximación y comprobar su validez deben de estar optimizados, de tal forma que sean lo más rápido posible. Este proyecto estará centrado en implementar y optimizar el Método del Eesiduo Equilibrado utilizando utilizando programa de cálculo computacional MATLAB. 1 Estructura del proyecto En este capítulo se expondrá la disposición de contenido que se tendrá a lo largo de este proyecto. En primer lugar, se hará una descripción sobre la base teórica del Método de los Elementos Finitos en la que de destacarán sus principales ideas características. Además se verá el problema que será considerado en todo el proyecto. Después, se explicará la base teórica del Método del Residuo Equilibrado, viendo cada una de las elecciones que se hacen para conseguir expresar la estimación del error como la superposición de la influencia de cada elemento. Con todo esto, se podrá generar un primer código de MATLAB que implemente cada uno de los pasos del Método del Residuo Equilibrado con el fin de ser una base sobre la que trabajar. Por último se estudiará el rendimiento del mismo y se irán realizando diferentes modificaciones que hasta llegar a un código definitivo, consecuencia del proceso de optimización. 3 2 Método de los Elementos Finitos En esta sección se hará una descripción del Método de los Elementos Finitos, que será el utilizado a la hora de calcular la solución aproximada, cuyo error se estimará con el Método del Residuo Equilibrado. 2.1 Introducción Como se nombró anteriormente, a la hora de resolver problemas se encuentran diversas ecuaciones algebráicas, intregrales, diferenciales o en derivadas parciales que no pueden ser resueltas de forma directa. Así, es necesario a recurrir a métodos numéricos para obtener una solución aproximada. A lo largo de la historia se han usado algunos como el Método de las Diferencias Finitas, con el que se discretiza el domino y obtiene una aproximación mediante la linealización de las ecuaciones diferenciales del problema en cada una de sus partes. De modo que se transforma una ecuación diferencial en otra algebraica. Otra opción utilizada por el Método de Galerkin y el de los Mínimos cuadrados fue la de obtener una solución como la combinación lineal de funciones, cjψj , estableciendo las funciones aproximantes, ψj, y siendo cjlas incognitas. [1] En base a esta última idea surge el Método de los Elementos Finitos, a mediados del siglo XX en la aplicación de análisis estructural de sistemas de barras. Con la particularidad de que, además de hacer una discretización del dominio, se hace la elección de las funciones aproximantes unas definidas a trozos, conocidas como funciones de pequeño soporte. Su uso se extrapoló a otros muchos campos como la Tranferencia de Calor, el Electromagnetismo o la Mecánica de Fluidos.[1] 2.2 Base teórica Inicialmente, el problema se define con la búsqueda de u(x), para cada cada xexistente en el dominio Ω, que satisface: D(u) = fpara x∈Ω u=uDpara x∈ΓD ∂ u ∂ n =gpara x∈ΓN            donde: •D(u) = f, ecuación diferencial que se tiene que cumplir en todos los puntos. 5 6Capítulo 2. Método de los Elementos Finitos •ΓD, frontera Dirichlet, en la que se conoce el valor de la incógnita. • Γ N , frontera Neumann, en la que el valor del flujo de u vendrá dada por la función g . •∂Ω, contorno del dominio. Se cumple que ∂Ω=ΓDSΓN. Para poder obtener una solución aproximada mediante el Método de los Elementos Finitos es necesario reformular el problema en formulación débil, es decir, se escriben de forma integral y las soluciones soluciones pertenecen a un espacio de Banach V , es decir, es un espacio normado y completo, que en resumidas cuentas se puede expresar como que cada vector contenido en el mismo tiene norma finita y que cualquier diferencia entre dos de sus vectores está en V. Quedando el problema expresado de la siguiente forma: u∈Vtal que B(u,v) = L(v),∀v∈V donde: •V, espacio de soluciones. •B ( u,v ), función que toma como entrada dos vectores de una misma longitud, con los que se obtiene otro con las mismas dimensiones. •L(v), función que transforma un vector en otro de su misma dimensión.[3] [2] El siguiente paso en el MEF es relizar la discretización del dominio, P, es decir, se dividirá el dominio Ωen un cierto número de partes, conocidas también como elementos, con formas definidas por líneas rectas, quedando unos vértices en cada uno de ellos que se conocen como nodos. Tras llevar a cabo este proceso se tiene un mallado, cuyas principales características son tener un número nTde elementos y nNnodos. Si se tiene un problema unidimensional, cada elemento será un segmento definido por los nodos de sus extremos. Mientras que en los problemas bidimensionales, los elementos suelen ser triángulos, con 3 nodos, o cuadriláteros, con 4 nodos. Figura 2.1 Ejemplo de una discretización con elementos triangulares. [2]. Una vez se tienen ésto, se puede definir la solución aproximada, ux, de igual forma que en el Método de Galerkin y el de los Mínimos cuadrados. La aproximación será la combinación lineal entre varias funciones, quedando definida por la elección de las funciones ψi, y los coeficientes ci, que están por determinar. De tal forma que: ux(x) = N X i=1 ciψi(x) La gran ventaja del MEF nace con la elección de estas funciones aproximantes psii, ya que, se toman funciones definidas a trozos asociadas a cada nodo, conocidas como funciones de pequeño soporte. De forma que son 1 en el nodo y 0 en en los otros nodos, es decir, únicamente serán no nulas en los elementos colindantes al nodo al que están 2.3 Problema considerado 7 asociadas. En la práctica, estas funciones son desde lineales hasta cúbicas. En consecuencia a esta elección, el valor de los coeficientes asociados a cada función será el valor de la solución aproximada en el nodo al que están asociadas. Las funciones de pequeño soporte se suelen ser simples, siendo lineales, cuadráticas o cúbicas. [1] Figura 2.2 Función de pequeño soporte lineal en caso unidimendional.[2]. Figura 2.3 Función de pequeño soporte lineal en caso bidimendional.[2]. Incluyendo esto a la formulación del problema se llega a encontrar Uique cumpla lo siguiente: nN X i=1 ciB(ψi,ψj) = l(ψj), j = 1,2...nN Es decir, se tiene un sistema de ecuaciones lineales de nNincógnitas con nNecuaciones, que puede ser representada de forma simplificada como K u =f. En la que Kse conoce como matriz de rigidez, en la Mecánica de Sólidos Continuos. Gracias a que las funciones son de pequeño soporte, el término que multiplica a cada ciserá no nulo cuando i=jo cuando i y j estén en un mismo elemento. Consecuentemente la matriz K será una matriz con una cantidad importante de elementos nulos, lo que hace que sea muy ventajoso definirla como una matriz dispersa. [2] Una vez se calculan esos coeficientes ya se tienen definidas todas las variables necesarias para conocer el valor de la solución aproximada en cada punto del dominio. 2.3 Problema considerado A la hora de trabajar en la optimización del funcionamiento del código del Método del Residuo Equilibrado se fijará un problema, que será resuelto mediante el MEF para varios mallados del dominio. El problema será bidimensional con una solución de una función escalar, por lo que 8Capítulo 2. Método de los Elementos Finitos se podría asimilar a un problema de transmisión de calor en una placa, en el que la incógnita es el valor de la temperatura. El dominio de este problema, Ω, será un cuadrado de lado de longitud 1 en el que se tiene la siguiente ecuación en derivadas parciales y condiciones de contorno. Ω = [0,1]×[0,1] −▽2u+u=fen Ω ∂ u ∂ n =gen ΓN u= 0 en ΓD siendo ▽2u=∂2u ∂x2+∂2u ž∂y2 f= 16x(1−x)+16y(1−y)+8x(1−x)y(1−y) g=(−8y(1−y)si x = 0 ox= 1 −8x(1−x)si y = 0 oy= 1 Este problema se solucionará con un programa de MATLAB dado, utilizando elementos triangulares y con funciones de pequeño soporte lineales. El mallado el número de elementos del mallado estará definido por el parámetro N, número de nodos que hay en cada lado de la frontera. Además se generarán de dos formas, una en la que todos los triángulos tienen el mismo tamaño, y otra en la que sus dimensiones son aleatorias, asimilándose más a las discretizaciones que se dan en la práctica, en la que el mallado es más fino en unas zonas que en otras. A continuación se expone un ejemplo de cada tipo, con su correspondiente solución. Figura 2.4 Mallado regular de 8 elementos. 2.3 Problema considerado 9 Figura 2.5 Mallado regular de 32 elementos. 16 Capítulo 3. Método del Residuo Equilibrado Figura 3.2 Ejemplo de nodo interior. [4]. Figura 3.3 Ejemplo de nodo fronterizo. [4]. y por los lados interiores µγ2 2,n +µγ2 1,n = 0 µγ3 3,n +µγ3 2,n = 0 . . . µγN N−1,n +µγN N−1,n = 0              Quedan por determinar 2 restricciones, una en el lado γ1y otra en el γN+1, que dependerán del tipo de frontera que sean los lados que llegan al nodo. Se pueden dar los siguientes tres casos: 2.a Neumann-Neumann. Los lados γ1 y γN+1 pertenecen a Γ N , por lo que se tienen las siguientes dos imposiciones para los momentos µγ1 1,n =Zγ1 gθnds µγN+1 N,n =Zγ1+N gθnds 2.b Neumann-Dirichlet. El lado γ1está en la frontera Neumann y γN+1 en la frontera Dirichlet, es válido para el caso de que sean al contrario, ya que sólo se tendrían que intercambiar las restricciones µγ1 1,n =Zγ1 gθnds µγN+1 N,n =sin restricciones 3.2 Descripción del método 17 2.c Dirichlet-Dirichlet.os lados γ1yγN+1 están contenidos en ΓD, entonces µγ1 1,n =sin restricciones µγN+1 N,n =sin restricciones Tabla 3.1 Tabla resumen de todos los casos y tipo de solución de cada uno.[4]. Incógnitas Restricciones Momentos Elementos Lados Lados Caso de flujo en equilibrio interiores exteriores Solución 12N N N 0No única 2.a 2N N N −1 2 Única 2.b 2N N N −1 1 Única 2.c 2N N N −1 0 No única Como se observa en la tabla, existe solución para todos los casos, aunque para dos de ellos la solución no es única. Es decir, se pueden dar infinitos casos, unos más cercanos en los que la diferencia entre los momentos reales y los aproximados, que son los que se calculan con los sistemas anteriores, puede llegar a variar mucho según la solución que se tome. Así que para conseguir que sean lo más parecidas posibles se tomará la que minimiza la diferencia entre ambas. Luego la función objetivo a minimizar será la siguiente 1 2 nT X K=1 X γ⊂∂K µγ K,n −˜µγ K,n2.(3.21) siendo µγ K,n los momentos reales y ˜µγ K,n los momentos aproximados. µγ K,n =Zγ θnnK·▽uds ˜µγ K,n =Zγ θnnK·▽uxds [4] Para encontrar los óptimos se utiliza la herramienta de los multiplicadores de Lagrange. Quedando la siguiente Lagrangiana L{˜µγ K,n},{λγ},{σK}=1 2 nT X K=1 X γ⊂∂K µγ K,n −˜µγ K,n2+ nT X K=1  ∆K(θn)−X γ⊂∂K µγ K,n + +X γ=∂K∩∂K′ λγ,n µγ K,n +µγ K′,n+X γ=∂K∩ΓN λγ,n µγ K,n −Zγ gθnds. Por comodidad se fija que λγ,n sea nulo en los lados de frontera Dirichlet debido a que no hay restricciones en los momentos de esos lados. De forma que a las restricciones de (3.19) se le suma µγ K,n −˜µγ K,n −σK,n +λγ,n = 0 (3.22) y λγ,n = 0 sobre γ⊂ΓD(3.23) 18 Capítulo 3. Método del Residuo Equilibrado Los multiplicadores de cada lado se puede expresar uniendo estas dos últimas restricciones junto a las de (3.19). Llegando a : λγ,n =           1 2σK,n +σK′,n + ˜µγ K,n + ˜µγ K′,nγ=∂K ∩∂K σK,n + ˜µγ K,n −Zγ gθnds γ =∂K ∩ΓN 0γ=∂K ∩ΓD Con estas relaciones, se puede expresar los momentos del flujo como función de σK,n y los momentos aproximados, que son conocidos al ser una entrada la solución aproximada de elementos finitos. µγ K,n =           1 2σK,n −σK′,n + ˜µγ K,n −˜µγ K′,nγ=∂K ∩∂K Zγ gθnds γ =∂K ∩ΓN σK,n + ˜µγ K,n γ=∂K ∩ΓD Para poder expresar todas las restricciones sobre los momentos del flujo según el multiplicador de Langrange se incluye esto último en la primera ecuación de (3.19). 1 2X γ=∂K∩∂K′σK,n −σK′,n+X γ⊂∂K∩ΓD σK,n =˜ ∆K(θn)∀K∈ Pn(3.24) donde ˜ ∆K(θn)=∆K(θn)−1 2X γ=∂K∩∂K′˜µγ K,n −˜µγ K′,n−X γ⊂∂K∩ΓD ˜µγ K,n −X γ⊂∂K∩ΓNZγ gθnds = =BK(ux,θn)−(f,θn)K−Z∂K ∂ux ∂nKθnds siendo ∂ux ∂nK=         1 2nK·{(▽ux)K+(▽ux)K′}sobre ∂K ∩∂K′ nK·(▽ux)Ksobre ∂K ∩ΓD gsobre ∂K ∩ΓN De esta forma, los multiplicadores de Langrange son las únicas incógnitas calcular para obtener los momentos, y consecuentemente el flujo. Éstos se obtienen de (3.24), llegando a los siguientes sistemas de ecuaciones en cada caso: 1. Nodo interior: 1 2         2−1... −1 −1 2 −1... 0 . . .. . . 0... −1 2 −1 −1... −1 2                 σ1 σ1 . . . σN−1 σN         =        ˜ ∆1(θn) ˜ ∆2(θn) . . . ˜ ∆N−1(θn) ˜ ∆N(θn)         2. Nodo fronterizo. Se dan los tres casos comentados anteriormente. 3.2 Descripción del método 19 2.a Neumann-Neumann: 1 2         1−1... 0 −1 2 −1... 0 . . .. . . 0... −1 2 −1 0... −1 1                 σ1 σ1 . . . σN−1 σN         =        ˜ ∆1(θn) ˜ ∆2(θn) . . . ˜ ∆N−1(θn) ˜ ∆N(θn)         2.b Neumann-Dirichlet: 1 2         3−1... 0 −1 2 −1... 0 . . .. . . 0... −1 2 −1 0... −1 1                 σ1 σ1 . . . σN−1 σN         =        ˜ ∆1(θn) ˜ ∆2(θn) . . . ˜ ∆N−1(θn) ˜ ∆N(θn)         2.c Dirichlet-Dirichlet: 1 2         3−1... 0 −1 2 −1... 0 . . .. . . 0... −1 2 −1 0... −1 3                 σ1 σ1 . . . σN−1 σN         =        ˜ ∆1(θn) ˜ ∆2(θn) . . . ˜ ∆N−1(θn) ˜ ∆N(θn)         [4] Se puede comprobar que las matrices de del vértice interior y Neumann-Neumann son singulares. Ambas tienen como espacio nulo el vector l= [1,1,...,1]T, así que se añade la restricción de que el vector σ sea ortogonal al espacio nulo, o lo que es lo mismo, que la suma de todas sus componentes sea nula. De esta forma sí se obtiene una solución única y se puede obtener el valor del flujo, con el que calcular el error. [4] Finalmente, aquí es donde se consiguen ver los puntos favorables de este método, ya que, en lugar de tener que lidiar con el problema completo para construir los flujos, se tienen unos sistemas definidos para cada nodo, que pueden llegar a tener matrices de coeficientes iguales, luego algunos de ellos se podrán llegar a resolver simultáneamente.[ 4 ] 4 Implementación en MATLAB Una vez se ha descrito el Método del Residuo Equilibrado, se procede a trazar su introducción en MATLAB. En primer lugar, se enumerarán cada una de las entradas que se obtienen desde el programa de elementos finitos. Después, se verán cada una de las partes del código, acompañadas de los razonamientos y expresiones que se usan en ellas, conformando un código sobre el que se optimizará más tarde. 4.1 Datos de entrada al código de MATLAB Para poder obtener una estimación del error, según este método, se reutilizan las siguientes matrices del programa de elementos finitos. Para poner expresar las dimensiones de las matrices se considera un mallado genérico de elementos triangulares compuesto por nTelementos y nNnodos. 4.1.1 Información del mallado Estas matrices definen completamente el mallado y proceden del programa de la obtención de la discretización del dominio, que se realiza justo antes de obtener la solución aproximada. •T, matriz de triangulación. Tiene una dimensión de 3×nTy contiene los índices de los nodos que forman cada elemento. Siendo la columna k la que contiene los índices de los nodos del elemento k,[n1,k,n2,K ,n3,k]T. T=   n1,1n1,2n1,nT n2,1n2,2... n2,nT n3,1n3,2n3,nT    •Z, matriz de coordenadas. Matriz de 2×nNen la se incorporan las coordenadas de todos los nodos, teniendo en la columna nla coordenada xeydel nodo n. Z="x1x2... xnN y1y2... ynN# •B, matriz de contorno. Matriz de 2×Lfront con la que se establecen los lados que componen la frontera indicando sus dos extremos, estando contenidos los índices de los nodos de los extremos del lado j de la frontera columna j . Su dimensión depende 21 22 Capítulo 4. Implementación en MATLAB del número de lados que compongan la frontera, cuyo número será la suma de los lados de frontera Dirichlet y los de frontera Neumann, Lfront =LN+LD. B=  nl,1nl,2... nl,nf ront nr,1nr,2... nr,nfront   •IN, índices de frontera Neumann. Vector de LNcomponentes con los índices de los lados de la frontera que son de tipo Neumann. IN = [γN,1,γN,2,...,γN,LN] •ID, índices de frontera Dirichlet. Vector de LDcomponentes con los índices de los lados de la frontera que son de tipo Dirichlet. ID = [γD,1,γD,2,...,γD,LD] •ite, matriz de vecindad. Tiene una dimensión de nT×3y se utiliza para cuál es el triángulo vecino a cada lado todos los elementos, o si es frontera, de la siguiente forma: ite(i,j) =    k > 0, el lado jdel triángulo ies común al triángulo k k < 0, el lado jdel triángulo ies el lado −kde la frontera •iteb , matriz complementaria de vecindad. Al igual que ite su tamaño es de nT× 3y ,como su nombre indica, añade información que no aporta la matriz de vecindad sobre cuál es el lado del triángulo vecino con el que se coincide. iteb(i,j) =    l > 0, el lado jdel triángulo ies el lado ldel triángulo k l= 0 , el lado jdel triángulo ies frontera, luego no hay vecino 4.1.2 Datos extraídos del cálculo de la solución aproximada Estas matrices se calculan en el programa de resolución del problema de elementos finitos, y se pueden reutilizar para obtener una estimación del error. Reducen el tiempo de la obtención del mismo, ya que no hay que volver a calcularlas •u, solución en los nodos. Vector de nNcomponentes que contiene el valor de la solución aproximada en cada nodo. u= [u1,u2,...,unN]T •Ae. Matriz de 3nT×ecompuesta de nTsubmatrices 3×3, una para cada elemento. Ae =      B1(θi,θj) B2(θi,θj) . . . BnT(θi,θj)       4.2 Procedimiento inicial en MATLAB 23 BK(θi,θj) =    BK(θ1,θ1)BK(θ1,θ2)BK(θ1,θ3) BK(θ2,θ1)BK(θ2,θ2)BK(θ2,θ3) BK(θ3,θ1)BK(θ3,θ2)BK(θ3,θ3)   •fe. Matriz de 3×nTque contiene las integrales sobre cada elemento de fθnpara todos los triángulos. fe =    (f,θ1)1(f,θ1)2(f,θ1)nT (f,θ2)1(f,θ2)2(f,θ2)nT (f,θ3)1(f,θ3)2(f,θ3)nT     •ge. Matriz de 2×LNque contiene la integral de gθnsobre todos los lados de frontera Neumann para sus dos extremos, siendo la columna j el valor de la integral para el lado de la frontera cuyo indice está en la posición jde IN. ge =     ZγN1 gθlds ZγN2 gθlds ... ZγNLN gθlds ZγN1 gθrds ZγN2 gθrds ... ZγNLN gθrds      donde los subíndices ryldenotan al nodo inicial y final de cada lado, respectivamente. 4.2 Procedimiento inicial en MATLAB Se describirán los pasos que se han seguido al programar en MATLAB con el fin de implementar el Método del Residuo Equilibrado. Este primer código, llamado versión 1.0, se elabora siguiendo la descripción del apartado anterior con la idea de poder entender cómo se introduce cada una de las variables del problema a MATALB, y el uso que se les da para obtener el resultado final. Todo ello sin tener en mente el tiempo que requieren para ser calculado, ya que, el código ideado aquí servirá como base, de estructura y de variables necesarias, para todas sus modificaciones posteriores. 1. Manipulación de los datos del mallado. Aquí se obtiene NT y Nn , número de triángulos y nodos, que serán necesarios para conocer las dimensiones de las matrices que se tendrán que crear más tarde. Estas constantes son fáciles de obtener mediante el uso del comando size sobre una de las filas de las matrices Tyz. Además se obtiene la longitud de cada lado de los triángulos y se guarda en Ck, una matriz de 3×nT, quedando en la columna jlas longitudes de los lados en este orden [ L1,L2,L3 ] T , tomando los lados de los triángulos como se ve en la Figura 4.1. Su cálculo es trivial, ya que, únicamente hay que calcular la longitud de un segmento. Cγ K=rxf−xi2+yf−yi2 siendo (xi,yi)y(xf,yf), las coordenadas los extremos del lado γdel triángulo K. 24 Capítulo 4. Implementación en MATLAB Figura 4.1 Triángulo Kgenérico. Código 4.1 Datos del Mallado. NT=length(T(1,:)); %Numero de triangulos Nn=length(z(1,:)); %Numero de nodos Ck=zeros(3,NT);%Longitud de los lados. Ck(L,k) lado L del triangulo k for k=1:NT %Para cada triangulo indices=T(:,k); %Indices de los nodos del triangulo k %- - Longitud de lados - - Ck(:,k)=[sqrt((z(1,indices(3))-z(1,indices(2)))^2... ...+(z(2,indices(3))-z(2,indices(2)))^2); sqrt((z(1,indices(1))-z(1,indices(3)))^2... ...+(z(2,indices(1))-z(2,indices(3)))^2); sqrt((z(1,indices(2))-z(1,indices(1)))^2... ...+(z(2,indices(2))-z(2,indices(1)))^2)]; end 2. Cálculo del gradiente de ux. Se guardará en la matriz gradux de 2×nT, con una columna para cada triángulo, siendo la componente de la primera fila la componente según Xy la otra fila la de Y. En primer lugar, se verá la expresión para obtener el valor del gradiente de la solución aproximada para un triángulo genérico. Para ello, se utiliza el siguiente cambio de variables para transformar todos los triángulos en uno rectángulo con vértices en (0,0),(1,0) y(1,1) en el nuevo sistemas de coordenadas. "x(ξ,η) y(ξ,η)#="x1 y1#+"x2−x1x3−x1 y2−y1y3−y1#"ξ η# La expresión de ux(ξ,η)sería ux(ξ,η)=[u1,u2,u3]    φ1(ξ,η) φ2(ξ,η) φ3(ξ,η)     4.2 Procedimiento inicial en MATLAB 25 siendo φnlas funciones base o de forma de cada nodo en el triángulo de vértices (0,0), (1,0) y(1,1). N(ξ,η) =     φ1(ξ,η) φ2(ξ,η) φ3(ξ,η)     =    1−ξ−η ξ η     Así que aplicando la regla de la cadena el gradiente sería de la siguiente forma:      ∂ux ∂x ∂ux ∂y      =     ∂ux ∂ξ ∂ξ ∂x +∂ux ∂η ∂η ∂x ∂ux ∂ξ ∂ξ ∂y +∂ux ∂η ∂η ∂y      = [u1,u2,u3]hNξξx+Nηηx,Nξξy+Nηηyi siendo Nξ= [−1,1,0]T,Nη= [−1,0,1]T ξx=y3−y1 Jηx=−y2−y1 J ξy=−x3−x1 Jηy=x2−x1 J J= (x2−x1)(y3−y1)−(x3−x1)(y2−y1) En definitiva, se tiene una expresión universal gracias a que el gradiente de esta forma. Así que se realizará un bucle con el que se obtendrán sus dos componentes para cada triángulo. Cabe destacar que el valor del gradiente es constante en cada elemento gracias a que las funciones de forma con las que se trabaja en este caso son lineales. Código 4.2 Gradiente de ux. gradux=zeros(2,NT);%Fila 1 u_x - - -Fila 2 u_y for k=1:NT %Para cada triangulo indices=T(:,k); %Indices de los nodos del triangulo k Nxi=[-1;1;0]; Neta=[-1;0;1]; denominador=(z(1,indices(2))-z(1,indices(1)))*(z(2,indices(3))-... ...z(2,indices(1)))-(z(1,indices(3))-z(1,indices(1)))*... ...(z(2,indices(2))- z(2,indices(1))); xix=(z(2,indices(3))-z(2,indices(1)))/denominador; etax=-1*(z(2,indices(2))-z(2,indices(1)))/denominador; xiy=-1*(z(1,indices(3))-z(1,indices(1)))/denominador; etay=(z(1,indices(2))-z(1,indices(1)))/denominador; gradux(1,k)= [u(indices(1)) u(indices(2)) u(indices(3))]*... ...(Nxi*xix+Neta*etax); gradux(2,k)= [u(indices(1)) u(indices(2)) u(indices(3))]*... ...(Nxi*xiy+Neta*etay); end 32 Capítulo 4. Implementación en MATLAB else %es frontera,hay que ver el tipo [~,frontera]=find(ID==frontera); if isempty(frontera) %tipo N-N sigmaden(1:N,n)=VerticeNeuNeu(deltagorron); elseif length(frontera)==1 %tipo N-D sigmaden(1:N,n)=VerticeDirNeu(deltagorron); else %tipo D-D sigmaden(1:N,n)=VerticeDirDir(deltagorron); end end end %---- %se asigna a sigma sigma=zeros(3,NT); for k=1:NT %cada triangulo indices=T(:,k); for n=1:3 %en cada nodo fila=find(Tsigma(:,indices(n))==k); sigma(n,k)=sigmaden(fila,indices(n)); end end function [ sigma ] = VerticeNeuNeu( TIndep ) %VERTICE DE FRONTERA NEUMANN-NEUMANN,singular entonces (N+1,1) ec extra del espacio nulo N=length(TIndep); if N==1 sigma=TIndep/1*2; else b=[TIndep; 0]; d0=[1 2*ones(1,N-2) 1]; d1=-1*ones(1,N-1); M=zeros(N+1,N); M(1:N,1:N)=diag(d0,0)+diag(d1,-1)+diag(d1,1); M(N+1,:)=ones(1,N); M=M/2; [L,U,P] = lu(M); %Ax=b PA=LU LUx=Pb Ly=Pb Ux=y y=L\(P*b); sigma=U\y; end end function [sigma]=VerticeInterior(TIndep) %VERTICE INTERIOR %singular entonces (N+1,1) ec extra del espacio nulo N=length(TIndep); 4.2 Procedimiento inicial en MATLAB 33 b=[TIndep; 0]; d0=2*ones(1,N); d1=-1*ones(1,N-1); Mint=zeros(N+1,N); Mint(1:N,1:N)=diag(d0,0)+diag(d1,-1)+diag(d1,1); Mint(1,N)=-1; Mint(N,1)=-1; Mint(N+1,:)=ones(1,N); Mint=Mint/2; [L,U,P] = lu(Mint); %Ax=b PA=LU LUx=Pb Ly=Pb Ux=y y=L\(P*b); sigma=U\y; end function [ sigma ] = VerticeDirDir( TIndep ) %VERTICE DE FRONTERA DIRICHLET-DIRICHLET %MATRIZ SIMETRICA DEFINIDA POSITIVA => CHOLESKY N=length(TIndep); if N==1 sigma=TIndep/3*2; else d0=[3 2*ones(1,N-2) 3]; d1=-1*ones(1,N-1); M=diag(d0,0)+diag(d1,-1)+diag(d1,1); M=M/2; R=chol(M); %A=R’R R’y=b Rx=y y=R’\TIndep; sigma=R\y; end end function [ sigma ] = VerticeDirNeu( TIndep ) %VERTICE DE FRONTERA DIRICHLET-NEUMANN %MATRIZ SIMETRICA DEFINIDA POSITIVA => CHOLESKY N=length(TIndep); if N==1 sigma=TIndep/3*2; else d0=[3 2*ones(1,N-2) 1]; d1=-1*ones(1,N-1); M=diag(d0,0)+diag(d1,-1)+diag(d1,1); M=M/2; R=chol(M); %A=R’R R’y=b Rx=y y=R’\TIndep; sigma=R\y; end 34 Capítulo 4. Implementación en MATLAB end 6. Construcción de los momentos de los flujos. Una vez se han calculado los valores de los multiplicadores y los momentos aproximados se puede tener el valor de µγ K,n, que será almacenado en la matriz mu de igual forma que se hace con muAprox, es decir, en una matriz (3×3nT)compuesta de nTsubmatrices de (3×3). Para explicar el procedimiento seguido se recuerda la expresión para obtener µγ K,n. µγ K,n =           1 2σK,n −σK′,n + ˜µγ K,n −˜µγ K′,nγ=∂K ∩∂K Zγ gθnds γ =∂K ∩ΓN σK,n + ˜µγ K,n γ=∂K ∩ΓD Al depender del lado y del elemento, se irá de lado a lado por cada uno de los elementos, o lo que es lo mismo, llenando una a una las columnas de la matriz que guarda el valor de los momentos del flujo. Se identificarán los lados de igual forma que en el caso del cálculo de ˜ ∆K,n y el triángulo vecino y su lado mediante ite eiteb. Debe de destacarse que en el caso de lado interior, las componentes de ˜µγ K′,n y las de σK′,n que se toman deben estar referidas a los nodos del triángulo K, es decir, se han de ordenar mediante la matriz Ordenacion, al igual que se hizo en el cálculo de los momentos del flujo. Mientras que cuando el lado sea frontera Neumann, se sigue el mismo procedimiento que en muAprox. Código 4.6 Valor del flujo. mu=zeros(3,3*NT); for k=1:NT for L=1:3 %ver tipo if ite(k,L)>0 %no es frontera Tvec=ite(k,L);%triangulo vecino Lvec=iteb(k,L);%lado del triangulo vecino sigmavec=sigma(Ordenacion(:,3*(L-1)+Lvec),Tvec); muAproxvec=muAprox(Ordenacion(:,3*(L-1)+Lvec),... ...3*(Tvec-1)+Lvec); mu(:,3*(k-1)+L)=0.5*(sigma(:,k)-sigmavec+... ..muAprox(:,3*(k-1)+L)-muAproxvec); elseif isempty(find(IN==-ite(k,L)))%es frontera, ver tipo %es frontera Dir mu(:,3*(k-1)+L)=sigma(:,k)+muAprox(:,3*(k-1)+L); else%es frontera Neu es dato LNeu=find(IN==-ite(k,L)); if L==1 mu(:,3*(k-1)+L)=[0;ge(:,LNeu)]; elseif L==2 mu(:,3*(k-1)+L)=[ge(2,LNeu);0;ge(1,LNeu)]; else mu(:,3*(k-1)+L)=[ge(:,LNeu);0]; end 4.2 Procedimiento inicial en MATLAB 35 end end end 7. Integral sobre los lados y estimación del error. Con todos los pasos anteriores ya se tienen todos los datos necesarios para definir la función del flujo sobre cada uno de los lados de los elementos, definida en la sección anterior como: gK|γ=2 hn2µγ K,l −µγ K,rθl+−µγ K,l +2µγ K,rθro El valor que interesa calcular sobre el flujo es su integral sobre los lados de todos los elementos, y teniendo en cuenta que es una función lineal, conociendo el valor del flujo en dos puntos y la longitud del lado se puede obtener el valor de la integral como la media entre los dos puntos multiplicado por la longitud del lado. Los flujos en los extremos se guardarán en la matriz gk de (2×3nT)con una columna para cada lado de cada elemento y una fila para el extremo inicial del lado y otra para el final. Se va llenando yendo de triángulo en triángulo y de lado en lado, definiendo gK en los extremos del lado de cada iteración, es decir, uno para cuando θl= 0 y el otro cuando θl= 0. En cuanto la integral se almacena en Intgk , un vector de nT componentes, cada una de ellas asociada a un elemento. Para obtener su valor se realiza con un bucle para barrer todos los triángulos y se va sumando la aportación de cada lado como la media del flujo de los extremos por la longitud de los lados. Código 4.7 Flujo y su integral sobre los lados. gk=zeros(2,3*NT); %fila 1 valor en nodo inical, fila 2 valor en nodo final for k=1:NT %lado 1 de 2 a 3 %lado 2 de 3 a 1 %lado 3 de 1 a 2 ThetaInicial=[2,3,1]; ThetaFinal=[3,1,2]; for L=1:3 gk(:,3*(k-1)+L)=2/Ck(L,k)*[2*mu(ThetaInicial(L),3*(k-1)+L)... ...-mu(ThetaFinal(L),3*(k-1)+L); -mu(ThetaInicial(L),... ...3*(k-1)+L)+2*mu(ThetaFinal(L),3*(k-1)+L)]; end end %----- %integral de gk en el contorno de todos los triangulos Intgk=zeros(1,NT); for k=1:NT for L=1:3 Intgk(k)=Intgk(k)+0.5*(gk(1,3*(k-1)+L)+... ...gk(2,3*(k-1)+L))*Ck(L,k); end end 5 Optimización del rendimiento 5.1 Rendimiento de la versión 1.0 En la sección anterior se ha planteado un primer código que implementa el Método del Residuo Equilibrado, en esta se estudiará los tiempos de cálculos y se intentarán reducir todo lo posible. En el desarrollo de esta sección se dividirán los tiempos de programa de MATLAB en los siguientes grupos: •Parte 1. Cálculo de ˜ ∆K,n. Incluye desde el paso 1 hasta el 4 de la sección anterior. •Parte 2. Obtención de σK,n, es en la que se resuelven todos los sistemas de ecuaciones que surgen en cada nodo para calculas los multiplicadores de Langrange. •Parte 3. Cálculo de la integral sobre el contorno de cada elemento del flujo. En ella se obtiene el valor de los momentos, con los que se reconstruye el flujo en cada lado. Estos tiempos se comprobarán para diferentes mallados, que se definirán por N, el número de divisiones de cada lado del rectángulo inicial, y la forma en la que se crea. Una de ellas será con un mallado en el que todos los elementos son triángulos rectángulos del mismo tamaño. Mientras que la otra será creada de forma aleatoria. Figura 5.1 Ejemplos de mallados. Como primera referencia de tiempos se tienen los de cálculo de la solución de elementos finitos para varios valores de N, que quedan recogidos a continuación en las Tablas 5.1 y 5.2, en las que se incluye también otros datos sobre la discretización del dominio en cuestión, como el número de elementos y de nodos. 37 38 Capítulo 5. Optimización del rendimiento En estas tablas, se observa que los tiempos del primer caso son más pequeños que los del segundo, esta diferencia se acentúa al crecer el refinamiento de la malla. Además, atendiendo al número de elementos de las mallas se ve que han tomado desde unas relativamente bastas, hasta unas con una cantidad de elementos varios órdenes de magnitud superior. Tabla 5.1 Tiempos de cálculo de la solución de elementos finitos para red regular. N20 80 160 240 330 450 nºde triángulos 800 12800 51200 115200 217800 405000 nºde nodos 441 6561 25921 58081 109561 203401 tsol (s) 0.4412 0.1751 0.8464 2.2335 4.2326 8.0890 Tabla 5.2 Tiempos de cálculo de la solución de elementos finitos para red generada aleatoriamente. N20 80 160 240 330 450 nºde triángulos 800 12800 51200 115200 217800 405000 nºde nodos 441 6561 25921 58081 109561 203401 tsol (s) 1.2973 0.6582 3.1478 8.1834 18.6325 27.9692 Con todo esto ya si se puede empezar a comparar tiempo con la primera versión del código que se tiene. Teniendo en cuenta que se reutilizan datos del problema de elementos finitos y el mallado es ya conocido, se espera tener unos tiempos mucho más pequeños que los de referencia. Entonces, se utiliza el código midiendo los tiempos de cada parte para poder hacer una comparación. Figura 5.2 Tiempos de cálculo con mallados regulares en V 1.0. 5.1 Rendimiento de la versión 1.0 39 Tabla 5.3 Tiempos de cálculo por partes para red regular en V 1.0 . N20 80 160 240 330 450 Parte 1(s) 0.0149 0.1020 0.5721 1.1652 2.1159 3.7622 Parte 2(s) 0.0152 0.6573 6.1568 22.9602 73.8507 311.7745 Parte 3(s) 0.0069 0.0923 0.4862 1.0292 1.7764 3.1908 Total(s) 0.0370 0.8516 7.2151 25.1546 77.7431 318.7275 Solución MEF(s) 0.4412 0.1751 0.8464 2.2335 4.2326 8.0890 Figura 5.3 Tiempos de cálculo con mallados irregulares en V 1.0. Tabla 5.4 Tiempos de cálculo por partes para red generada aleatoriamente. N20 80 160 240 330 450 Parte 1(s) 0.0147 0.1462 0.5504 1.3347 2.3189 3.1528 Parte 2(s) 0.0301 1.0479 10.2125 44.0317 156.0442 397.1398 Parte 3(s) 0.0089 0.1271 0.4672 1.0494 2.1089 2.6992 Total(s) 0.0537 1.3212 11.2301 46.4158 160.4720 402.9918 Solución MEF(s) 1.2973 0.6582 3.1478 8.1834 18.6325 27.9692 En ambas situaciones se tienen situaciones parecidas, los tiempos de la estimación del error es bastante superior al del cálculo de la solución aproximada. Ésto significa que el primer código, no es para nada válido para ser utilizado para mallados con muchos elementos, ya que, su rendimiento en ambos es muy bajo y tiene un gran margen de mejora a priori. Tener tiempos de cálculo más altos que los de elementos finitos es inaceptable, al estar reutilizando muchos datos . La aportación más grande al tiempo de cálculo, especialmente para mallados con muchos elementos, viene de la parte 2, en la que se resuelven los sistemas de ecuaciones. Según se ha planteado inicialmente, se resuelven nNsistemas de ecuaciones, de modo que cuanto mayor sea la Nmucho mayor es el número de nodos, y consecuentemente mayor es el número de matrices de coeficientes que hay que crear y de sistemas que se resuelven. 40 Capítulo 5. Optimización del rendimiento 5.2 Versión 1.1 Esta versión nace con la siguiente pregunta, ¿ se pueden reducir el número de matrices de coeficientes que se crean? La respuesta es sí, debido a que siempre son iguales para los nodos que están en el mismo número de elementos y que son del mismo tipo, como se vio en la sección del Método del Residuo Equilibrado. AX1=b1 AX2=b2 . . . AXn=bn ≡A  X1X2... Xn  =  b1b2... bn   Luego un buen cambio en el planteamiento para la nueva versión, llamada versión 1.1, es ir resolviendo a la vez todos los sistemas de los nodos que están en un mismo número de triángulos. Siendo solamente necesario crear cuatro matrices de coeficientes, una para cada tipo de nodo, para cada número de triángulos en los que está un nodo, cifra que determina la dimensión de las matrices. El número de triángulos en los que esta cada nodo se puede conocer siguiendo el mismo procedimiento que se hacía para encontrar su valor máximo, mediante el comando sparse, quedando almacenado en el vector NTri siendo la componente i el número de triángulos en los que está el nodo de índice i . Además su valor máximo es el valor NTrimax, que determinará el tamaño de sigmaden, al igual que en la versión anterior. NTri=sparse(T,ones(3*NT,1),ones(3*NT,1)); Introduciendo este idea al obtener la matriz sigma, se pasa de crear nNveces las matrices de coeficientes y resolver su sistemas a 4NTrimax, considerando que siempre hay al menos un nodo de cada tipo para cada número de triángulos en los que está un nodo. Generalmente, el valor máximo de NTri no tiene por qué ser mucho mayor en redes finas que en mallados poco refinados, sin embargo el número de nodos si es mucho mayor, luego este cambio debería aportar una reducción en el tiempo de cálculo para redes con un número cuantioso de nodos. Con este planteamiento, al calcular σK,n hay que hacer un bucle desde 1 hasta NT rimax, en el que se identifica los nodos que están el número de triángulos en cuestión mediante el comando find. Luego de esto, hay que separarlos por tipos. Para ello, se hace un bucle en por cada uno de ellos, identificándose de la misma forma que se hacía en la versión inicial, llenando unos vectores para cada tipo en los que se tienen el índice del nodo y unas matrices para los triángulos en los que están y el nodo que son cada uno dentro de ese triángulo. Además de tener en una matriz de la misma estructura con el valor de ˜ ∆K,n correspondiente. Después de ésto, se crean las matrices de coeficientes de cada tipo y se resuelven para todos los nodos iguales a la vez, aprovechando que la matriz de coeficientes es igual. La matriz de coeficientes de un tipo se crea siempre que haya al menos un nodo del mismo. Es decir, la solución de cada tipo será una matriz compuesta por columnas que contienen la solución de cada nodo. Por último, se asignan los valores de σa su correspondiente posición en la matriz en la que se guarda σK,n, de la misma forma que se hacía en la versión 1.0 . Código 5.1 Calculo de los multiplicadores de Lagrange. 5.2 Versión 1.1 41 %Calculo sigma sigma=zeros(size(T)); NTri=sparse(T,ones(3*NT,1),ones(3*NT,1)); nmax=max(NTri);% for j=1:nmax Nodosj=find(NTri==j); %nodos que estan en j triangulos NodosInt=zeros(length(Nodosj),1); %nodos interiores nInt=0;%numero de nodos TInt=zeros(j,length(Nodosj));% triangulos nkInt=TInt; %nodos de los triangulos deltaInt=TInt; %termino independiente NodosNeuNeu=NodosInt; %nodos neu-neu nNeuNeu=0; TNeuNeu=TInt; nkNeuNeu=TInt; deltaNeuNeu=TInt; NodosNeuDir=NodosInt;%nodos neu-dir nNeuDir=0; TNeuDir=TInt; nkNeuDir=TInt; deltaNeuDir=TInt; NodosDirDir=NodosInt; %nodos dir-dir nDirDir=0; TDirDir=TInt; nkDirDir=TInt; deltaDirDir=TInt; %tipo de nodo for n=Nodosj’ [filas,columnas]=find(T==n); deltagorron=zeros(j,1); for k=1:j%vector con los indices de los triangulos deltagorron(k)=deltagorro(filas(k),columnas(k)); end %es frontera? [~,frontera]=find(B==n); if isempty(frontera) %es interior nInt=nInt+1; NodosInt(nInt)=n; TInt(:,nInt)=columnas; nkInt(:,nInt)=filas; 48 Capítulo 5. Optimización del rendimiento denominador=((X(2,:)-X(1,:)).*(Y(3,:)-Y(1,:))-... ...(X(3,:)-X(1,:)).*(Y(2,:)-Y(1,:))); xix=(Y(3,:)-Y(1,:))./denominador; etax=-(Y(2,:)-Y(1,:))./denominador; xiy=-(X(3,:)-X(1,:))./denominador; etay=(X(2,:)-X(1,:))./denominador; Además se usa el comando kron, que realiza el producto de Kronecker a los matrices dadas. Esta operación, sobre una matriz A ( m×n )y otra B ( p×q ), obtiene una matriz de (mp×nq)con la siguiente relación: A⊗B=   a11B ... a1nB . . ..... . . am1B ... amnB    donde ⊗denota producto de Kronecker. [5] Con el producto de Kronecker sobre Nξ y ξx , siendo ξT x un vector con un valor columna con una componente para cada elemento, se tiene una matriz de (3 ×nT )de la siguiente forma: Nξ⊗ξT x=   −ξx(1) −ξx(2) ... −ξx(nT) ξx(1) ξx(2) ... ξx(nT) 0 0 ... 0   Luego si se hace un sumatorio de todas sus filas, se llega a un vector de nT componentes que sería equivalente al producto que se hacía en cada bucle. Este razonamiento se extrapola para los otros términos con los que se calcula el gradiente. Llegando a las siguientes líneas de código para calcular el gradiente sin usar ningún bucle. Código 5.6 Obtención del gradiente. gradux=zeros(2,NT); gradux(1,:)=sum(U.*(kron(Nxi,xix(:)’)+kron(Neta,etax(:)’))); gradux(2,:)=sum(U.*(kron(Nxi,xiy(:)’)+kron(Neta,etay(:)’))); 3. Obtención del momento del flujo aproximado, ˜µγ K,n. Este paso se realizaba con dos bucles, uno para cada triángulo y otro para cada lado de los mismos, luego conseguir eliminarlos sería muy ventajoso. Para conseguir hacerlo , en primer lugar se guarda en una matriz (3×nT)el producto escalar entre la normal exterior a cada lado y el gradiente de la solución aproximada en ese elemento mediante al producto elemento a elemento entre ambos. Siendo la componente ( i,j )de esa matriz el producto escalar en el elemento jentre su gradiente y la normal exterior al lado (i,j). Para calcular el flujo se utiliza el comando repmat con el que se pueden crear matrices a raíz de repetir las filas y columnas de una matriz o un vector. Utilizando esto sobre la matriz de las longitudes, expresada como un vector fila, la matriz de los productos escalares, de igual modo que la anterior, y repitiendo nT veces la matriz IntThetan , se 5.3 Versión 1.2 49 llega a poder calcular de la misma forma de la que se hacía lado a lado. Veáse el ejemplo para un triángulo K, tomando las submatrices correspondientes. ˜µγ K,n =ZLj nj K·▽ux|Kθids = =      n1 K·▽ux|Kn2 K·▽ux|Kn3 K·▽ux|K n1 K·▽ux|Kn2 K·▽ux|Kn3 K·▽ux|K n1 K·▽ux|Kn2 K·▽ux|Kn3 K·▽ux|K       .∗      C1|KC2|KC3|K C1|KC2|KC3|K C1|KC2|KC3|K       .∗      01 2 1 2 1 201 2 1 2 1 20       donde .∗ denota el producto elemento a elemento. En la práctica se realiza ese cálculo a la vez para todos los triángulos, ya que las matrices que se utilizan son de (3 × 3 nT ), es decir, están compuestas por nTsubmatrices con esa forma. Código 5.7 Cálculo del momento del flujo aproximado. IntThetan=[0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0 ]’; %En eta xi gradun=gradux(1,:).*nx+gradux(2,:).*ny; %esta bien muAprox=repmat(Ck(:)’,3,1).*repmat(IntThetan,1,NT).*... ...repmat(gradun(:)’,3,1); 4. Calculo de ˜ ∆K,n. Recordando la ecuación de este término: ˜ ∆K(θn) = BK(ux,θn)−(f,θn)K−Z∂K ∂ux ∂nKθnds Al calcular sus dos primeros términos se usaba un bucle para ir elemento a elemento, que se elimina fácilmente en el caso de ( f,θn ) K , ya que, en el fondo sólo se creaba una matriz igual a fe, que es dato del problema. Mientras que el caso de BK(ux,θn)tiene algo de complejidad. Originalmente se hacia el producto de la submatriz de Ae ,correspondiente a cada elemento, con el vector con el valor de uxen sus nodos. Para poder hacer esta multiplicación matricial simultáneamente en todos los triángulos se sabe que esta operación se compone de multiplicar elemento a elemento cada fila de Ae|Kcon ux|Ky después sumarlos, de forma que quedan asignados en la fila correspondiente de Ae|K. Luego si se consigue hacer el primer paso con una multiplicación elemento a elemento sólo habría que realizar un sumatorio de filas. Con este razonamiento se crea una matriz (3nT×3) compuesta por nTmatrices que tienen 3 filas repetidas, iguales a ux|K, mediante el uso de kron sobre U, matriz con valores nodales de ux, y un vector de unos. De esta forma, se tiene una matriz con la que multiplicar elemento a elemento con Ae . Después quedaría hacer un sumatorio de las filas y se tendría un vector de 3nTcomponentes, compuesto por nTvectores que contienen [BK(ux,θ1),BK(ux,θ2),BK(ux,θ3)]T. Véase a continuación el ejemplo de ésto en triángulo K.       BK(ux,θ1) BK(ux,θ2) BK(ux,θ3)       =      Ae11|KAe12|KAe13|K Ae21|KAe22|KAe33|K Ae31|KAe32|KAe33|K       .∗      u1u2u3 u1u2u3 u1u2u3       50 Capítulo 5. Optimización del rendimiento Este vector se transforma en una matriz Bk como la que se tenía en versiones anteriores mediante el comando reshape, que va llenando una matriz de dimensiones dadas con las componentes de un vector dado. fek=fe; Bk=reshape(sum(Ae.*kron(U’,ones(3,1)),2),3,NT); Ahora sólo queda el último término, que como se vio anteriormente depende del tipo de cala lado. Antes se calculada con dos bucles y dentro de ellos se utilizaban varios if, por lo que eliminarlos podría reducir los tiempos de cálculo cuando hay un grande número de elementos. Aquí también se consiguen eliminar todos los bucles, con la particularidad de que se tienen que realizar 3 cálculos diferentes, uno para cada tipo de lado. Es decir, para todos los lados del mismo tipo se calcula simultáneamente, siendo necesario obtener la información para definir todos los lados cada tipo. La matriz TiposLados se ha creado para eso, mediante la búsqueda del valor asignado a cada tipo en esa matriz con find, se obtiene I , un vector con los triángulos en los que está cada lado, y J , vector con el lado que es en cada triángulo que se encuentra en ese lugar de I. En el caso del lado interior, hay que realizar una semisuma entre el vector de momentos aproximados del flujo de cada lado con el opuesto de su triángulo vecino en cada uno de ellos. Se hará primero la parte del triángulo en cuestión y después se le añadirá la aportación del vecino. En primer lugar habrá que buscar en TiposLados el valor 0, consiguiendo los correspondientes IyJ. Con éstos se crea una matriz auxiliar compuesta por los momentos aproximados, que son vectores columna de 3 componentes, referentes a los lados que definen los vectores IyJ. Se utiliza sparse asignar estos valores en su posición correspondiente de Intuxnk, matriz (3×nT)en la que se que almacena los valores de la integral para cada triángulo. La razón para elegir este comando es la misma que cuando se generó la matriz NTri, porque al querer sumar todos los lados de cada triángulo en una misma columna, si se da un valor a una componente que ya tiene, éstas se suman. Para asignar los valores almacenados en la matriz auxiliar multiplicados por 0.5, sólo habrá que crear la matriz de filas del mismo tamaño, compuesta por la repetición de estos vectores [1,2,3]T, que indican a las funciones base que se han tomado para calcular cada momento; y otra matriz con las columnas, que será la repetición 3 veces del vector I. Código 5.8 Primera parte de la integral del producto escalar en los lados internos. [I,J]=find(TiposLados==0); aux=muAprox(:,(I-1)*3+J); fila=repmat((1:3)’,1,length(aux(1,:))); columna=repmat(I’,3,1); Intuxnk=0.5*sparse(fila,columna,aux,3,NT); En estos tipos de lados todavía queda incluir la aportación de los triángulos vecinos, éstos y sus correspondientes lados quedarán definidos por los valores de IyJen la matriz ite eiteb. Se extrae los vectores Ivec yJvec extrayendo esos datos de las matrices citadas. Con éstos datos se repite el proceso hecho anteriormente, con la particularidad que el signo de la matriz auxiliar será negativo y de que hay que ordenar el orden de los índices de los nodos, para que coincidan en los dos triángulos colindantes, mediante el uso de la matriz Ordenacion. Es decir, la matriz que determina las filas en las que se asignan los valores de la matriz auxiliar serán columnas 5.3 Versión 1.2 51 de Ordenacion, elegidas en función de cuál sea el valor del lado y lado vecino, almacenados en JyJvec respectivamente. Con todo ello se suma la nueva matriz dispersa a la que se tenía de antes. Código 5.9 Segunda parte de la integral del producto escalar en los lados internos. %L1 de k %L2 de k %L3 de k Ordenacion=[ 1 2 3 2 3 1 3 1 2; %theta 2 de k 312 123 231;%theta 2 de k 231 312 123];%theta 3 de k %L1,L2,L3 de k Ivec=ite((J-1)*NT+I); Jvec=iteb((J-1)*NT+I); aux=muAprox(:,(Ivec-1)*3+Jvec); %hay que cambiar el orden fila=Ordenacion(:,3*(J-1)+Jvec); Intuxnk=Intuxnk-0.5*sparse(fila,columna,aux,3,NT); Para el caso de los Dirichlet, sería prácticamente igual que la primera parte del caso anterior, con la diferencia de que se busca el valor 1en la matriz TiposLados y que la matriz auxiliar no se multiplica por 0.5. Código 5.10 Integral del producto escalar en los lados Dirichlet. [I,J]=find(TiposLados==1); aux=muAprox(:,(I-1)*3+J); fila=repmat((1:3)’,1,length(aux(1,:))); columna=repmat(I’,3,1); Intuxnk=Intuxnk+sparse(fila,columna,aux,3,NT); Sólo queda añadir la aportación de los lados de frontera Neumann. Se busca en la matriz TiposLados el valor 2, la matriz auxiliar será la matriz ge a la que se le añade una fila de ceros. Para poder asignar los valores a sus correspondientes funciones nodales y elemento hace falta que el orden de de los vectores IyJcoherentes a como están dispuestos los datos en ge. Primero se introduce la matriz OrdenLNeu, cuya función es la equivalente a la de Ordenacion, ordenando las componentes de ge en función del índice de cada lado. OrdenLNeu =   231 312 123   En la matriz ge los datos están dispuestos por columnas desde el lado de frontera Neumann con menor indice hasta el de mayor. Para conseguir que IyJestén dispuestos de esa forma, se extraen estos índices en el vector LNeu de la matriz ite los índices de lado de frontera asignados a cada valor de esta pareja de vectores. Después se obtiene el vector que los ordena mediante el comando sort sobre LNeu . Finalmente se repite el proceso del caso anterior, estando I y J re-ordenados. Y finalmente, se obtiene el valor de ˜ ∆K,n sumando sus tres términos. Código 5.11 Integral del producto escalar en los lados Neumann. 52 Capítulo 5. Optimización del rendimiento [I,J]=find(TiposLados==2); ge=[ge;zeros(1,length(ge(1,:)))]; OrdenLNeu=[2 3 1; 3 1 2; %Componente de ge 1 2 3]; %L1 L2 L3 LNeu=-ite((J-1)*NT+I); [~,orden]=sort(LNeu); aux=ge; fila=OrdenLNeu(:,J(orden)); columna=repmat(I(orden)’,3,1); Intuxnk=full(Intuxnk+sparse(fila,columna,aux,3,NT)); deltagorro=Bk-fek-Intuxnk; 5. Obtención de los multiplicadores de Langrange σK,n. Aquí es donde se gasta el mayor tiempo de cálculo para mallados finos, como se vio anteriormente, así que reducir las operaciones que realiza MATLAB resulta muy beneficioso para los tiempos. Al igual que con los lados, se calculará para cada tipo de nodos, siendo este el único paso en el que no se pueden eliminar todos los bucles. Ésto se debe a que realizar un bucle para crear las matrices de coeficientes correspondientes es imprescindible. Lo primero será crear un vector que contiene los índices de los nodos para cada tipo mediante el uso de find sobre la matriz TiposNodos , buscando su valor correspondiente. Antes de comenzar a ir resolviendo por cada tipo de nodo, se crean unas matrices con una información complementaria a NTri en cuya fila j se guardan los triángulos en los que está el nodo jy su índice en éstos. Ambas son matrices dispersas con posiciones vacías, luego para obtener un vector referente a un nodo con componentes no nulas se puede utilizar el comando nonzeros. Código 5.12 Clasificación de los nodos. NodosInt=find(TiposNodos(:,1)==0)’; NodosNeuNeu=find(TiposNodos(:,1)==1)’; NodosNeuDir=find(TiposNodos(:,1)==2)’; NodosDirDir=find(TiposNodos(:,1)==3)’; sigma=zeros(3,NT); NTri=sparse(T,ones(3*NT,1),ones(3*NT,1)); nmax=max(NTri);% Numero mas alto de triangulos en los que puede estar un nodo %triangulos en los que esta cada nodo TNodos=sparse(T,repmat(1:NT,3,1),repmat(1:NT,3,1)); %indices en los triangulos NNodos=sparse(T,repmat(1:NT,3,1),repmat((1:3)’,1,NT)); Con esto ya sí se puede comenzar a resolver los sistemas. Se empieza con los nodos interiores. Se hace un bucle que determina el número de triángulos en los que va a estar los nodos que se tomen, en este caso se empieza por un valor de 3 porque es imposible 5.3 Versión 1.2 53 que un nodo interior esté en un número menor. Geométricamente, todos los ángulos internos de un triángulo debe de sumar 180º, entonces al menos serán necesarios 3 para llevar a completar 360º. Para cada vuelta en el bucle, se buscan los nodos interiores que están en el número j de triángulos y sus índices se guardan en el vector Nodosj . Si hay al menos un nodo, se sigue el mismo proceso que en la versión anterior, pero con otros órdenes al obtenerse todos los datos de los nodos al mismo tiempo. Ahora, se extraen los valores de los triángulos en los que está cada nodo y sus índices en ellos directamente de las correspondientes filas de las matrices TNodos y NNodos , a lo que se le suma el uso de nonzeros. Siendo el resultado en cada matriz es un vector de que se contiene de varios subvectores de j componentes que contienen la información de cada nodo. Entonces para tenerlos de una forma igual a la de versiones anteriores se vuelve a utilizar reshape. Con los términos independientes pasa exactamente lo mismo. El siguiente paso es construir la matriz correspondiente a este tipo de nodo y a jy resolverlo , obteniéndose la matriz sigmaj . Finalmente, se introducen los valores en su posición correcta de sigma igual que en la versión 1.1. Código 5.13 Cálculo de los multiplicadores de Lagrange en nodos interiores. for j=3:nmax %nodos interiores que estan en j triangulos Nodosj=find(NTri(NodosInt)==j); %nodos que estan en j triangulos nj=length(Nodosj); if nj>0 Nodosj=NodosInt(Nodosj);%indices Tj=reshape(nonzeros(TNodos(Nodosj,:)’),j,length(Nodosj)); nkj=reshape(nonzeros(NNodos(Nodosj,:)’),j,length(Nodosj)); deltaj=reshape(deltagorro(nkj+3*(Tj-1)),j,length(Nodosj)); d0=2*ones(1,j); d1=-1*ones(1,j-1); M=zeros(j+1,j); M(1:j,1:j)=diag(d0,0)+diag(d1,-1)+diag(d1,1); M(1,j)=-1; M(j,1)=-1; M(j+1,:)=ones(1,j); M=M/2; b=[deltaj(:,1:nj); zeros(1,nj)]; sigmaj=M\b; sigma(nkj+(Tj-1)*3)=sigmaj; end end Para los otros tipos de nodos sería un proceso similar, cambiando el vector de nodos sobre el que se busca, la matriz de coeficientes y haciendo que j empiece desde 1, dado que ese caso se da en los nodos que se encuentran en las esquinas del dominio. 54 Capítulo 5. Optimización del rendimiento Código 5.14 Cálculo de los multiplicadores de Lagrange en nodos fronterizos. for j=1:nmax Nodosj=find(NTri(NodosNeuNeu)==j); nj=length(Nodosj); if nj>0 Nodosj=NodosNeuNeu(Nodosj); Tj=reshape(nonzeros(TNodos(Nodosj,:)’),j,length(Nodosj)); nkj=reshape(nonzeros(NNodos(Nodosj,:)’),j,length(Nodosj)); deltaj=reshape(deltagorro(nkj+3*(Tj-1)),j,length(Nodosj)); d0=[1 2*ones(1,j-2) 1]; d1=-1*ones(1,j-1); M=zeros(j+1,j); if j==1 M=1/2; sigmaj=M\deltaj; else M(1:j,1:j)=diag(d0,0)+diag(d1,-1)+diag(d1,1); M(j+1,:)=ones(1,j); M=M/2; b=[deltaj(:,1:nj); zeros(1,nj)]; sigmaj=M\b; end sigma(nkj+(Tj-1)*3)=sigmaj; end end for j=1:nmax Nodosj=find(NTri(NodosNeuDir)==j); nj=length(Nodosj); if nj>0 Nodosj=NodosNeuDir(Nodosj); Tj=reshape(nonzeros(TNodos(Nodosj,:)’),j,length(Nodosj)); nkj=reshape(nonzeros(NNodos(Nodosj,:)’),j,length(Nodosj)); deltaj=reshape(deltagorro(nkj+3*(Tj-1)),j,length(Nodosj)); d0=[3 2*ones(1,j-2) 1]; d1=-1*ones(1,j-1); 5.3 Versión 1.2 55 if j==1 M=3/2; else M=diag(d0,0)+diag(d1,-1)+diag(d1,1); M=M/2; end b=deltaj(:,1:nj); sigmaj=M\b; sigma(nkj+(Tj-1)*3)=sigmaj; end end for j=1:nmax Nodosj=find(NTri(NodosDirDir)==j); nj=length(Nodosj); if nj>0 Nodosj=NodosDirDir(Nodosj); Tj=reshape(nonzeros(TNodos(Nodosj,:)’),j,length(Nodosj)); nkj=reshape(nonzeros(NNodos(Nodosj,:)’),j,length(Nodosj)); deltaj=reshape(deltagorro(nkj+3*(Tj-1)),j,length(Nodosj)); if j==1 M=3/2; else d0=[3 2*ones(1,j-2) 3]; d1=-1*ones(1,j-1); M=diag(d0,0)+diag(d1,-1)+diag(d1,1); M=M/2; end b=deltaj(:,1:nj); sigmaj=M\b; sigma(nkj+(Tj-1)*3)=sigmaj; end end 6. Construcción de los momentos de los flujos, µγ K,n. Siguiendo un procedimiento equivalente al del cálculo de los flujos aproximados usando los vectores IyJde cada caso, añadiendo la influencia de los multiplicadores de Langrange en los lados interiores y Dirichlet. Código 5.15 Cálculo de los momentos del flujo. 56 Capítulo 5. Optimización del rendimiento mu=zeros(3,3*NT); [I,J]=find(TiposLados==0); aux=sigma(:,I)+muAprox(:,3*(I-1)+J); fila=repmat((1:3)’,1,length(aux(1,:))); columna=repmat((3*(I-1)+J)’,3,1); mu(fila+3*(columna-1))=0.5*aux; Ivec=ite((J-1)*NT+I); Jvec=iteb((J-1)*NT+I); aux=sigma(:,Ivec)+muAprox(:,3*(Ivec-1)+Jvec); fila=Ordenacion(:,3*(J-1)+Jvec); mu(fila+3*(columna-1))=mu(fila+3*(columna-1))-0.5*aux; [I,J]=find(TiposLados==1); aux=sigma(:,I)+muAprox(:,3*(I-1)+J); fila=repmat((1:3)’,1,length(aux(1,:))); columna=repmat((3*(I-1)+J)’,3,1); mu(fila+3*(columna-1))=aux; LNeu=-ite((J-1)*NT+I); [~,orden]=sort(LNeu); aux=ge; I=I(orden); J=J(orden); fila=OrdenLNeu(:,J); columna=repmat((3*(I-1)+J)’,3,1); mu(fila+3*(columna-1))=aux; 7. Integral del flujo sobre los lados. Cuando se reconstruyen los valores del flujo para los extremos del flujo se consigue eliminar los bucles para cada triángulo y lado. Aquí se crean matrices repetidas con los números de los nodos iniciales y finales de cada lado, igual que el caso anterior, pero remetidas nT veces para que se pueda calcular en todos los triángulos a la vez. Al igual que con el cálculo de ˜µγ K,n , se repite una matriz con las longitudes de los lados, con las que dividir elemento a elemento en este caso. De forma que se realiza este cálculo para todos los triángulos con una misma operación:   gl|K g2|K = 2  2µ1 K,2−µ1 K,32µ2 K,3−µ2 K,12µ3 K,1−µ3 K,2 −µ1 K,2+2µ1 K,3−µ2 K,3+2µ2 K,1−µ3 K,1+2µ3 K,2 ./  C1|KC2|KC3|K C1|KC2|KC3|K  donde ./ denota la división elemento a elemento. Código 5.16 Reconstrucción del flujo. ThetaInicial=repmat([2,3,1],1,NT); ThetaFinal=repmat([3,1,2],1,NT); 5.3 Versión 1.2 57 gk=2./repmat(Ck(:)’,2,1).*[2*mu(3*(0:3*NT-1)+ThetaInicial)-... ...mu(3*(0:3*NT-1)+ThetaFinal); -mu(3*(0:3*NT-1)+ThetaInicial)+2*mu(3*(0:3*NT-1)+ThetaFinal)]; Por último se calculan todas las integrales sobre todos los lados. Recordando que gk es una matriz de (2 × 3 nT ), en la que las columnas 1 , 4 , 7 , 10 ... están referidas a los lados 1 de todos los triángulos, las 2,5,8,11... a los lados 2 y las 3,6,9,12... a los lados 3. Entonces se puede calcular la integral en cada uno de ellos como la mitad de la suma del producto elemento a elemento del vector con la media del flujo para los lados 1 con el vector con las longitudes de los lados 1, que sería la primera fila de Ck . Quedando el valor de la integral en todo el contorno del elemento como la suma de todos los lados . Código 5.17 Integral del flujo del flujo. Intgk=0.5*(sum(gk(:,1:3:3*NT-2)).*Ck(1,:)+... ...sum(gk(:,2:3:3*NT-1)).*Ck(2,:)+sum(gk(:,3:3:3*NT)).*Ck(3,:)); Como se ha visto se han se ha pasado de tener varios bucles en cada parte a únicamente tener 3, barriendo desde 1 o 3 hasta valores que generalmente serán a lo sumo del orden de 15 o 20, mucho menor que la gran cantidad de bucles que corrían decenas de miles de veces con mellados muy refinados. Para ver si la mejora que se obtiene es importante se volverán a medir los tiempos de cálculo y se compararán entre las distintas versiones. Figura 5.6 Comparación de tiempos con mallados regulares. Índice de Tablas 3.1 Tabla resumen de todos los casos y tipo de solución de cada uno.[4] 17 5.1 Tiempos de cálculo de la solución de elementos finitos para red regular 38 5.2 Tiempos de cálculo de la solución de elementos finitos para red generada aleatoriamente 38 5.3 Tiempos de cálculo por partes para red regular en V 1.0 39 5.4 Tiempos de cálculo por partes para red generada aleatoriamente 39 5.5 Comparación entre los tiempos de cálculo de la parte 2 para red regular 45 5.6 Comparación entre los tiempos de cálculo de la parte 2 para red generada aleatoriamente 45 5.7 Comparación entre los tiempos totales de cálculo para redes regulares 58 5.8 Comparación entre los tiempos totales de cálculo para redes generadas aleatoriamente 58 5.9 Tiempos por partes de la Versión 1.2 para red regular 59 5.10 Tiempos por partes de la Versión 1.2 para redes irregulares 60 65 Índice de Códigos 4.1 Datos del Mallado 24 4.2 Gradiente de ux25 4.3 Cálculo del momento del flujo aproximado 26 4.4 Términos independientes de los sistemas 29 4.5 Obtención de los multiplicadores de Lagrange 31 4.6 Valor del flujo 34 4.7 Flujo y su integral sobre los lados 35 [ 40 5.1 Calculo de los multiplicadores de Lagrange 40 5.2 Datos del Mallado 46 5.3 Matriz de tipos de nodo 47 5.4 Matriz de tipos de lados 47 5.5 Cambio de variables 47 5.6 Obtención del gradiente 48 5.7 Cálculo del momento del flujo aproximado 49 5.8 Primera parte de la integral del producto escalar en los lados internos 50 5.9 Segunda parte de la integral del producto escalar en los lados internos 51 5.10 Integral del producto escalar en los lados Dirichlet 51 5.11 Integral del producto escalar en los lados Neumann 51 5.12 Clasificación de los nodos 52 5.13 Cálculo de los multiplicadores de Lagrange en nodos interiores 53 5.14 Cálculo de los multiplicadores de Lagrange en nodos fronterizos 53 5.15 Cálculo de los momentos del flujo 55 5.16 Reconstrucción del flujo 56 5.17 Integral del flujo del flujo 57 67 Bibliografía [1] Federico París, "Capítulo 12: El Método de los Elementos Finitos", Teoría de la Elasticidad, 1996. [2] Endre Süli, Lecture Notes on Finite Element Methods for Partial Differential Equations, 2019. [3] WIKIPEDIA , Método de los Elementos Finitos. [Consulta: 20 de julio 2021]. [Consultarla]. [4] Mark Ainsworth, J.Tinsley Oden, A posteriori error estimation in finite element analysis, 1997. [5] WIKIPEDIA , Producto de Kronecker. [Consulta: 26 de julio 2021]. [Consultarla]. 69