Estructuras de datos tipo grafo en los algoritmos de refinamiento y desrefinamiento basados en la bisección por el lado mayor. Aplicaciones

UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA DEPARTAMENTO DE MATEMÁTICAS TESIS DOCTORAL JOSÉ PABLO SUÁREZ RIVERO Las Palmas de Gran Canaria, 2001 ESTRUCTURAS DE DATOS TIPO GRAFO EN LOS ALGORITMOS DE REFINAMIENTO Y DESREFINAMIENTO BASADOS EN LA BISECCIÓN POR EL LADO MAYOR. APLICACIONES Título de la tesis: ESTRUCTURAS DE DATOS TIPO GRAFO EN LOS ALGORITMOS DE REFINAMIENTO Y DESREFINAMIENTO BASADOS EN LA BISECCIÓN POR EL LADO MAYOR. APLICACIONES Thesis title: GRAPHBASED DATA STRUTURES FOR REFINEMENT AND DEREFINEMENT ALGORITHMS BASED ON THE LONGEST EDGE BISECTION. APPLICATIONS UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA Departamento de Matem´aticas Departamento de Cartograf´ıa y Expresi´on Gr´afica en la Ingenier´ıa TESIS DOCTORAL ESTRUCTURAS DE DATOS TIPO GRAFO EN LOS ALGORITMOS DE REFINAMIENTO Y DESREFINAMIENTO BASADOS EN LA BISECCI ´ ON POR EL LADO MAYOR. APLICACIONES. JOS´ EPABLO SU´ AREZ RIVERO OCTUBRE DE 2001 Tesis Doctoral: ESTRUCTURAS DE DATOS TIPO GRAFO EN LOS ALGORITMOS DE REFINAMIENTO Y DESREFINAMIENTO BASADOS EN LA BISECCI ´ ON POR EL LADO MAYOR. APLICACIONES. Programa de Doctorado: Simulaci´on Num´erica en Ciencia y Tecnolog´ıa Bienio: 1997-1999 Departamento: Departamento de Matem´aticas Universidad: Universidad de Las Palmas de Gran Canaria Las Palmas de Gran Canaria, a de de 2001 Autor Director Jose Pablo Su´arez Rivero Dr. ´ Angel Plaza de La Hoz Agradecimientos Quisiera agradecer a mi director de tesis Dr. ´ Angel Plaza por la dedicaci´on e inter´es que ha mostrado desde el primer momento en toda mi actividad investigadora. Sus acertadas sugerencias y confianza depositada en m´ı han hecho posible convertir en realidad un proyecto de tesis como el realizado. La sinton´ıa en nuestro trabajo de investigaci´on sigue en la actualidad. Asimismo, al profesor Dr. Graham F. Carey quien en mi ´ultima estancia en el TICAM (Texas Institute for Computational and Applied Mathematics) de la Universidad de Texas, mostr´o un gran inter´es en mi trabajo de investigaci´on que lo ha fortalecido. Los densos meses de trabajo bajo su direcci´on han sido claves para el t´ermino de la tesis. A la profesora Mar´ıa-Cecilia Rivara de la Universidad de Chile quien tambi´en siempre ha estado colaborando y aportando con su experiencia y saber numerosas sugerencias e ideas contempladas en esta tesis. Resumen Los algoritmos de generaci´on de mallas han sido objeto de numerosas investigaciones en los ´ultimos a˜nos. Estos algoritmos constituyen herramientas b´asicas en los m´etodos num´ericos, como en el m´etodo de elementos finitos o en la aproximaci´on lineal a trozos de funciones, as´ı como en la industria gr´afica, como gr´aficos por ordenador, dise˜no y modelado de s´olidos y superficies, visualizaci´on, CAD/CAM, etc. Cada vez m´as se utilizan en entornos de Simulaci´on Num´erica en Ciencia y Tecnolog´ıa, en donde es cr´ıtico que los algoritmos sean eficientes en t´erminos de tiempo y espacio de c´omputo. En los algoritmos de refinamiento y desrefinamiento de mallas es crucial disponer de una buena estructura de datos que haga eficiente la ejecuci´on de los mismos. Ello se debe a las complejas operaciones que surgen en la partici´on de elementos, lo cual se acent´ua cuando la dimensi´on del problema aumenta y la cantididad de elementos crece r´apidamente. Esta tesis investiga una clase de estructuras de datos basadas en grafos que permite dise˜nar algoritmos eficientes de refinamiento y desrefinamiento de mallas y propone algunas aplicaciones. El trabajo est´a estructurado como sigue: en el cap´ıtulo primero se hace una introducci´on general de la generaci´on de mallas, incluyendo ideas, conceptos, algoritmos y t´ecnicas y aplicaciones. Asimismo, se presentan algunas ideas b´asicas de dise˜no de estructuras de datos. Se concluye el cap´ıtulo se˜nalando los objetivos propuestos en esta tesis. El cap´ıtulo segundo se dedica al estudio de las particiones geom´etricas de elementos simpliciales en 2D y 3D y de los algoritmos de refinamiento y desrefinamiento m´as importantes de la literatura, especialmente los basados en el esqueleto. i ii El cap´ıtulo tres hace una revisi´on de las estructuras de datos geom´etricas y espaciales m´as importantes en la representaci´on de mallas. El cap´ıtulo cuatro presenta una nueva estructura de datos basada en grafos que de forma natural se ajusta a los algoritmos de refinamiento y desrefinamiento basados en el esqueleto de Plaza y Carey. Se presentan las propiedades matem´aticas y computacionales de dicha estructura de datos. Asimismo, se proporcionan versiones nuevas de los algoritmos de refinamiento en dimensi´on dos y tres y del algoritmo de desrefinamiento en dimensi´on dos. Finalmente, se formulan y demuestran diversas propiedades de la partici´on en cuatro tri´angulos por el lado mayor. En el cap´ıtulo quinto se presentan aplicaciones: en la secci´on primera, modelos digitales del terreno y generalizaci´on de terrenos; en la segunda secci´on, niveles de detalle en gr´aficos por ordenador y V RML; y finalmente, en la secci´on tercera se incorpora el algoritmo de refinamiento en 2D a un c´odigo de elementos finitos comercial y se resuelve un problema el´ıptico no lineal. El cap´ıtulo seis se destina a presentar las conclusiones resultantes de esta tesis y se proponen diversas l´ıneas de trabajo futuro. Finalmente, se presenta un anexo que recoge las tres publicaciones m´as importantes que han dado lugar las investigaciones llevadas a cabo en la tesis. Abstract In the last years, mesh generation algorithms have been largely investigated. These algorithms are basic tools in numerical methods, as in the finite element method or in the pieceway linear approximation to a given function. The algorithms are also relevant in the graphic industry, as in computer graphics, surfaces and solids modeling and design, visualization, CAD/CAM, etc. Furthermore, in an increasing way, the algorithms are getting used in Science and Technology Numerical Simulation; in this wide topic it is specially critical to provide algorithms with efficient perfomance in terms of time and space. In the refinement and coarsening algorithms, an important point is to have a suitable data structure which makes it possible that the algorithms perform in an effcicient way. This is due to the complexity of operations in the elements partitions which turn in a critical issue when the dimension increases and the number of elements grows rapidly. The present thesis investigates a new class of data structures which are based on the graph theory. The data structures are used to design new efficient versions of refinement and coarsening algorithms and applications. The outline of the treatment is as follows: chapter one introduces mesh generation, including basic ideas, concepts, algorithms and techniques, and applications. Some preliminaries about data structures are given and the goals of the thesis close the chapter. Second chapter studies the geometric partitions for simplicial elements in both 2D and 3D. It is also discussed the refinement and coarsening algorithms which are more relevant in the literature, specially the algorithms based on the skeleton. iii ´ INDICE DE FIGURAS xi 3.12. Esquema de acceso indexado . . . . . . . . . . . . . . . . . . . . 95 3.13. Esquema de acceso Hash . . . . . . . . . . . . . . . . . . . . . . 96 3.14. ´ Arbol Cuaternario . . . . . . . . . . . . . . . . . . . . . . . . . 99 3.15. ´ Arbol Kd de dos claves . . . . . . . . . . . . . . . . . . . . . . . 101 4.1. Grafo de incidencia para un tetraedro . . . . . . . . . . . . . . . 105 4.2. Representaci´on gr´afica del esqueleto . . . . . . . . . . . . . . . . 106 4.3. Configuraciones posibles para el grafo del (n−1)-esqueleto . . . 109 4.4. Grafo del esqueleto y vector T. . . . . . . . . . . . . . . . . . . 110 4.5. Grafo del 1-esqueleto en 2D . . . . . . . . . . . . . . . . . . . . 112 4.6. Existencia de bucles en el grafo del 1-esqueleto . . . . . . . . . . 112 4.7. Patrones del G1para los tres tipos de tetraedros . . . . . . . . . 114 4.8. Configuraciones posibles de G1para los tres tipos de tetraedros 115 4.9. 2-esqueleto para un tetraedro . . . . . . . . . . . . . . . . . . . 116 4.10. Componentes con nodo com´un de G2 1yG2 2. . . . . . . . . . . . 118 4.11. Niveles de refinamiento 1 y 2 del ejemplo test . . . . . . . . . . 124 4.12. Nivel de refinamiento 3 del ejemplo test . . . . . . . . . . . . . . 125 4.13. Ejemplo de refinamiento . . . . . . . . . . . . . . . . . . . . . . 129 4.14. La condici´on de desrefinamiento . . . . . . . . . . . . . . . . . . 131 4.15. Mallas 1 y 2 del test del algoritmo 2D-SBD . . . . . . . . . . . . 134 4.16. Mallas 3 y 4 del test del algoritmo 2D-SBD . . . . . . . . . . . . 135 4.17. (a) Malla equilibrada, (b) malla semi-equilibrada . . . . . . . . . 137 4.18. Extensi´on por conformidad y LEEP de tri´angulos vecinos . . . . 139 4.19. 4T-LE en un tri´angulo acut´angulo . . . . . . . . . . . . . . . . . 141 4.20. (a) 4T-LE en un tri´angulo obtus´angulo (b) Partici´on 4T-LE . . 142 4.21. (a) Tri´angulos terminales estables (b) Partici´on 4T-LE . . . . . 146 4.22. Mallas de los problemas test . . . . . . . . . . . . . . . . . . . . 148 4.23. Refinamiento global 4T-LE y evoluci´on de M1 y M2. Malla Delaunay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 4.24. Refinamiento global 4T-LE y evoluci´on de M1 y M2. Malla pentagonal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 4.25. Grado de equilibrio. Malla pentagonal y Malla Delaunay . . . . 153 4.26. Malla inicial tipo Delaunay y nivel de refinamiento uno. Refinamiento global 4T-LE . . . . . . . . . . . . . . . . . . . . . . 154 4.27. Malla Delaunay. Niveles de refinamiento dos y tres. Refinamiento global 4T-LE . . . . . . . . . . . . . . . . . . . . . . . . . . 155 xii ´ INDICE DE FIGURAS 4.28. Malla pentagonal. Refinamiento global 4T-LE. Malla inicial y nivel de refinamiento uno . . . . . . . . . . . . . . . . . . . . . . 156 4.29. Malla pentagonal. Refinamiento global 4T-LE. Nivel de refinamiento dos y cuatro . . . . . . . . . . . . . . . . . . . . . . . 157 5.1. Representaci´on inicial del terreno de Agaete . . . . . . . . . . . 165 5.2. Malla 2D y curvas de nivel del terreno de Agaete . . . . . . . . 166 5.3. Representaci´on con error de 1 metro del terreno de Agaete . . . 167 5.4. Malla 2D y curvas de nivel con error de 1 metro del terreno de Agaete . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 5.5. Representaci´on con error de 6 metros del terreno de Agaete . . . 169 5.6. Malla 2D y curvas de nivel con error de 6 metros del terreno de Agaete . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 5.7. Representaci´on con error de 15 metros del terreno de Agaete . . 171 5.8. Malla 2D y curvas de nivel con error de 15 metros del terreno de Agaete . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 5.9. Representaci´on con error de 36 metros del terreno de Agaete . . 173 5.10. Malla 2D y curvas de nivel con error de 36 metros del terreno de Agaete . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 5.11. Gr´afica del error en metros frente al n´umero de nodos . . . . . . 175 5.12. Gr´afica del SNR en decibelios frente al n´umero de nodos . . . . 175 5.13. Tiempo de obtenci´on de las mallas frente al nivel de error . . . . 175 5.14. Representaci´on inicial del terreno de G´aldar . . . . . . . . . . . 176 5.15. Representaci´on con error de 10 metros del terreno de G´aldar . . 177 5.16. Malla 2D y curvas de nivel con error de 10 metros del terreno de G´aldar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178 5.17. Representaci´on con error de 15 metros del terreno de G´aldar . . 179 5.18. Malla 2D y curvas de nivel con error de 15 metros del terreno de G´aldar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 5.19. Representaci´on con error de 40 metros del terreno de G´aldar . . 181 5.20. Malla 2D y curvas de nivel con error de 40 metros del terreno de G´aldar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 5.21. Representaci´on con error de 90 metros del terreno de G´aldar . . 183 5.22. Malla 2D y curvas de nivel con error de 90 metros del terreno de G´aldar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 5.23. Gr´afica del error en metros frente al n´umero de nodos . . . . . . 185 ´ INDICE DE FIGURAS xiii 5.24. Gr´afica del SNR en decibelios frente al n´umero de nodos . . . . 185 5.25. Tiempo de obtenci´on de las mallas frente al nivel de error . . . . 185 5.26. Puntos de vistas y niveles de detalle . . . . . . . . . . . . . . . . 190 5.27. Malla 2D y representaci´on VRML con nivel de detalle 1 . . . . . 191 5.28. Malla 2D y representaci´on VRML con nivel de detalle 2 ( 10 % error) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192 5.29. Malla 2D y representaci´on VRML con nivel de detalle 3 ( 30 % error) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193 5.30. Malla 2D y representaci´on VRML con nivel de detalle 4 ( 60 % error) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194 5.31. Diagrama en bloques de PDE Toolbox . . . . . . . . . . . . . . 196 5.32. Soluci´on inicial de la ecuaci´on el´ıptica no lineal . . . . . . . . . 200 5.33. Soluci´on en el nivel de refinamiento 1 . . . . . . . . . . . . . . . 201 5.34. Soluci´on en el nivel de refinamiento 2 . . . . . . . . . . . . . . . 202 5.35. Soluci´on en el nivel de refinamiento 3 . . . . . . . . . . . . . . . 203 5.36. Soluci´on en el nivel de refinamiento 4 . . . . . . . . . . . . . . . 204 ´ Indice de cuadros 2.1. Prueba num´erica para la triangulaci´on similar 4T . . . . . . . . 72 2.2. Medidas estad´ısticas I . . . . . . . . . . . . . . . . . . . . . . . 74 2.3. Medidas estad´ısticas II . . . . . . . . . . . . . . . . . . . . . . . 74 3.1. Dispositivos de almacenamiento . . . . . . . . . . . . . . . . . . 94 4.1. Datos de la prueba 1 . . . . . . . . . . . . . . . . . . . . . . . . 123 4.2. Secuencia de tri´angulos distintos obtenidos por la aplicaci´on iterativa de la partici´on 4T-LE . . . . . . . . . . . . . . . . . . . 144 4.3. Malla Delaunay. M´etricas M1 y M2 . . . . . . . . . . . . . . . . 151 4.4. Malla pentagonal. M´etricas M1 y M2 . . . . . . . . . . . . . . . 152 4.5. Malla Delaunay. Grado de equilibrio . . . . . . . . . . . . . . . . 153 4.6. Malla Pentagonal. Grado de equilibio . . . . . . . . . . . . . . . 153 4.7. Malla Delaunay. Medida de ´angulos . . . . . . . . . . . . . . . . 158 4.8. Malla Pentagonal. Medida de ´angulos . . . . . . . . . . . . . . . 158 5.1. Niveles de detalle para el terreno de G´aldar . . . . . . . . . . . 190 xv Cap´ıtulo 1 Introducci´on y Objetivos La generaci´on de mallas es el proceso de discretizaci´on de un dominio f´ısico en subdominios m´as peque˜nos (elementos). Las mallas se encuentran en diversas ´areas del an´alisis num´erico dentro de la Ciencia e Ingenier´ıa aplicada. Ejemplos de ´estas son: ajuste de curvas y superficies en la teor´ıa de aproximaci´on, integraci´on num´erica, diferencias finitas, vol´umenes finitos y elementos finitos. Todas tienen en com´un la soluci´on num´erica de ecuaciones diferenciales en derivadas parciales, cuyo proceso en general puede ser descrito de la siguiente forma: 1. Definici´on inicial del dominio proporcionando algunos puntos. 2. Generaci´on de la malla y etiquetado de nodos. 3. Soluci´on num´erica del problema en el dominio discretizado. 4. Refinamiento o redistribuci´on de la malla y salto al paso 3 para mejorar la soluci´on. 5. Presentaci´on de la malla y soluci´on (visualizaci´on). No obstante, la generaci´on de mallas se ha convertido actualmente en un tema interdisciplinar que aparece en muchas otras ´areas. Citamos por ejemplo la industria gr´afica, y m´as concretamente, gr´aficos por ordenador, dise˜no y modelado de s´olidos y superficies, visualizaci´on, CAD/CAM, etc. Asimismo la geometr´ıa y la topolog´ıa computacional hacen obligada visita al tema de la generaci´on de mallas. Otras ´areas m´as aplicadas como cartograf´ıa por ordenador y sistemas de informaci´on geogr´aficos (GIS) han venido a recurrir a 1 2Cap´ıtulo 1. Introducci´on y Objetivos t´ecnicas de mallado y afines para la representaci´on y visualizaci´on de datos espaciales. Fundamentalmente, los dominios 2D pueden ser divididos en tri´angulos o cuadril´ateros mientras que los dominios 3D en tetraedros o hexaedros. Idealmente la forma y distribuci´on de elementos en el dominio son realizados por los algoritmos de generaci´on de mallas, autom´aticos o no autom´aticos. Ser´an autom´aticos, cuando la generaci´on de mallas s´olo requiere de unos pocos par´ametros iniciales y a continuaci´on se realiza el proceso, no requiriendo la participaci´on del usuario m´as que al inicio. Ser´an no autom´aticos en caso contrario, es decir, el proceso de la generaci´on de la malla requiere de decisiones y par´ametros externos. En las d´ecadas recientes, el m´etodo de elementos finitos ha llegado a ser un soporte fundamental en el dise˜no y an´alisis de aplicaciones industriales. Cada vez m´as, dise˜nos complejos y de gran magnitud est´an siendo simulados usando el m´etodo de elementos finitos, lo que ha hecho aumentar la popularidad e inter´es por los algoritmos de mallado autom´atico. El problema de la generaci´on autom´atica de mallas considera la mejor descripci´on de un dominio geom´etrico, sujeto a un conjunto dado de nodos iniciales, especificaci´on de tama˜nos para los elementos y criterios de la forma de los elementos. La geometr´ıa del dominio se compone habitualmente de v´ertices, curvas, superficies y s´olidos, que en el mejor de los casos se describen por paquetes de CAD o de modelado de s´olidos. Existen dos tipos fundamentales de mallas: estructuradas y no estructuradas. En las mallas estructuradas la conexi´on entre elementos obedece a un patr´on fijo; en las no estructuradas no existe tal patr´on, el n´umero de nodos vecinos de uno dado es variable. Las mallas obtenidas por una transformaci´on de coordenadas, de un cuadrado en dimensi´on dos ´o un cubo en dimensi´on tres, son ejemplos de mallas estructuradas. En ellas el nodo de ´ındices (i, j, k) siempre tiene como vecino por la izquierda al (i−1, j, k) y por la derecha al (i+ 1, j, k), etc. 1.1. Algoritmos de generaci´on de mallas 3 Las mallas estructuradas tienen ciertas ventajas sobre las no estructuradas. Son m´as simples y tambi´en m´as convenientes para su uso en el m´etodo de diferencias finitas. Requieren menos memoria para su almacenamiento en ordenador y es m´as f´acil controlar la forma y tama˜no de los elementos. La gran desventaja que ofrecen las mallas estructuradas es su poca flexibilidad para adaptarse a dominios con geometr´ıas complicadas. Se han desarrollado numerosas t´ecnicas para encontrar las transformaciones de coordenadas adecuadas: transformaci´on conforme, m´etodos algebraicos, etc. Incluso con estas t´ecnicas es imposible encontrar una transformaci´on que se adapte satisfactoriamente a la forma de un dominio complicado. Existen multitud de problemas formulados en dominios cuya forma no es adecuada para una discretizaci´on mediante mallas estructuradas. En estos casos es necesario utilizar mallas no estructuradas. La generaci´on autom´atica de mallas no estructuradas es un campo relativamente reciente y que en poco tiempo ha producido grandes avances. Esta disciplina utiliza conceptos de la geometr´ıa computacional, y por consiguiente del ´algebra lineal y de la teor´ıa de poliedros entre otros campos, y sus aplicaciones se encuentran principalmente en el m´etodo de los elementos finitos. Es frecuente que, los elementos que componen las mallas no estructuradas sean s´ımplices, esto es, tri´angulos en dimensi´on dos y tetraedros en dimensi´on tres. La utilizaci´on de s´ımplices es debida, entre otras cosas, a su capacidad para llenar el espacio y para adaptarse a los contornos del dominio sin que se produzcan discontinuidades. La malla formada por s´ımplices recibe el nombre de triangulaci´on y en dimensi´on tres teselaci´on, o tambi´en triangulaci´on en dimensi´on tres. Es sabido que es posible descomponer un dominio poli´edrico (poligonal en dimensi´on dos) en tetraedros (tri´angulos), aunque esta descomposici´on no es ´unica y puede requerir la inserci´on de puntos adicionales. 1.1. Algoritmos de generaci´on de mallas Presentamos un breve resumen de los principales tipos de algoritmos utilizados para la generaci´on de mallas no estructuradas. La generaci´on de mallas no estructuradas a partir de tri´angulos y tetraedros es una de las t´ecnicas m´as 4Cap´ıtulo 1. Introducci´on y Objetivos utilizadas actualmente. La mayor´ıa de las t´ecnicas que se utilizan actualmente se pueden englobar en las tres principales categor´ıas: 1. Quadtree-octree 2. Delaunay 3. Avance Frontal Aunque hay ciertamente un grado de complejidad mayor cuando trabajamos en 3D que cuando lo hacemos en 2D, los algoritmos que aqu´ı se exponen en gran parte son aplicables para la generaci´on de mallas ya trabajemos en dos o tres dimensiones. 1.1.1. Quadtree-octree La t´ecnica Octree fue desarrollada inicialmente por el grupo de Mark S. Shephard en el Polythechnic Institute de Rensselaer, Nueva York [102, 103, 92] durante la d´ecada de los ochenta, y es una extensi´on del caso bidimensional conocido como Quadtree. Figura 1.1: Estructura de datos tipo quadtree En el m´etodo Octree, el concepto de quadtree (cuadrados que se dividen en cuatro, resultando una estructura de datos arb´orea, como se pone de manifiesto en la figura 1.1) es reemplazado por hexaedros -o cubosque se dividen entre cero y ocho hijos. Ahora el proceso de divisi´on es m´as complejo que en dimensi´on dos. En concreto para obtener la conformidad de la malla resultante, 1.1. Algoritmos de generaci´on de mallas 5 todos los elementos se dividen finalmente en tetraedros. Sea Ω un dominio acotado de IR3con frontera ∂Ω poligonal. Los pasos del algoritmo, de la misma naturaleza que en dimensi´on dos, son los siguientes: Paso 1: Se genera un paralelep´ıpedo inicial, que contiene todos los puntos de la frontera del dominio. El m´etodo consiste en dividir, partiendo de la malla inicial, cada caja o hexaedro-padre en algunos otros hexaedros hijos (entre cero y ocho) mediante un proceso recursivo hasta que, se obtiene la malla final. Esta malla final cumple el criterio de que cada caja final contiene como mucho un v´ertice de la frontera del dominio Ω. Este criterio permite controlar el tama˜no de cada elemento. Paso 2: La partici´on resultante del paso anterior se puede suavizar definiendo los hexaedros hijos de tal forma que, como mucho, haya un punto intermedio en cada cara de los hexaedros. Paso 3: Se realiza un an´alisis de las caras cuadrangulares teniendo en cuenta: Cara cuadrangular externa que no contenga ning´un punto de la frontera de Ω es eliminada. Cada caja interna que no contenga ning´un punto produce un hexaedro que es dividido en tetraedros, con el fin de obtener la conformidad de la malla. Para cada caja que contenga un 00trozo00 de frontera, se calculan los puntos intersecci´on, y el tetraedro es dividido. Se define en este punto el contorno final de la malla. Sin entrar en mucho detalle, merece la pena se˜nalar que la creaci´on de los tetraedros, partiendo de la malla de hexaedros se realiza en dos pasos: Cada cara cuadrangular es dividida en tri´angulos como se muestra en la figura 1.2 (b), siguiendo ciertos 00patrones00 semejantes al caso bidimensional, de forma que la superficie formada por las caras cuadrangulares (esqueleto) se hace conforme. Se dividen internamente los hexaedros en tetraedros de forma congruente con la divisi´on realizada previamente en las caras (ver figura 1.2 (c)). 6Cap´ıtulo 1. Introducci´on y Objetivos (a) Detalle de una malla de hexaedros (b) Divisi´on de las caras cuadrangulares en tri´angulos (c) Configuraci´on final Figura 1.2: Divisi´on de un cubo por octree 1.1. Algoritmos de generaci´on de mallas 13 Figura 1.10: Triangulaci´on de Delaunay, no ´unica P Figura 1.11: Inserci´on de un nuevo punto 14 Cap´ıtulo 1. Introducci´on y Objetivos tri´angulos forman un pol´ıgono denominado pol´ıgono de inserci´on, como se ve en la figura 1.12. Figura 1.12: Formaci´on de un pol´ıgono de inserci´on La nueva triangulaci´on de Delaunay estar´a formada por todos los tri´angulos de la anterior triangulaci´on, exteriores al pol´ıgono de inserci´on, m´as los tri´angulos formados al unir dicho punto con los v´ertices del pol´ıgono de inserci´on. P Figura 1.13: Nueva triangulaci´on de Delaunay La triangulaci´on de Delaunay re´une tambi´en algunas caracter´ısticas que la hacen adecuada para la aplicaci´on del M´etodo de Elementos Finitos (M.E.F.). Adem´as de la capacidad de refinamiento, la triangulaci´on de Delaunay posee propiedades geom´etricas ´optimas: de todas las triangulaciones de un conjunto de puntos del plano, la de Delaunay es la que hace m´aximo el m´ınimo ´angulo de cualquier tri´angulo. Sin embargo, no existe una generalizaci´on de 1.1. Algoritmos de generaci´on de mallas 15 esta propiedad para dimensi´on mayor que dos, al menos en t´erminos de alguna medida angular de los s´ımplices. En dimensi´on mayor o igual a dos hay propiedades que hacen referencia al tama˜no de las esferas que contienen a los s´ımplices de la triangulaci´on, pero que no nos dan una idea tan clara de la calidad de la triangulaci´on como la que da la propiedad anterior [24]. A pesar de ello, los resultados experimentales revelan que las mallas tridimensionales, construidas a partir de la triangulaci´on de Delaunay, dan buenas medidas de calidad cuando la distribuci´on de puntos es suficientemente regular. Una de las principales dificultades que plantea la utilizaci´on de la triangulaci´on de Delaunay para generar mallas es conseguir que la triangulaci´on respete la frontera del dominio de definici´on, esto es, que no haya ning´un tetraedro que corte su superficie. Cuando esto ocurre decimos que la triangulaci´on es conforme con la frontera del dominio. En dimensi´on dos se resuelve este problema, bien imponiendo que la triangulaci´on contenga todas las aristas que definen el contorno del dominio (Constrained Delaunay Triangulation), o bien, colocando puntos sobre las aristas del dominio en posiciones adecuadas de manera que ´estas sean la uni´on de aristas de la triangulaci´on (Conforming Delaunay Triangulation). Tambi´en en dimensi´on tres hay algoritmos (ver por ejemplo [31, 100]) basados en el intercambio de aristas y caras que consiguen la conformidad con la frontera del dominio. Una vez que se ha definido una malla conforme con la frontera del dominio, una buena estrategia para adaptar la malla, asegurando que en todo momento se mantiene dicha conformidad, ser´ıa utilizar m´etodos de refinamiento/desrefinamiento de mallas encajadas, ´o bien utilizar la triangulaci´on de Delaunay y el refinamiento mediante la t´ecnica de insertar los nuevos v´ertices en los puntos medios de las aristas [84]. Por otra parte, la realizaci´on pr´actica de algoritmos que generen la triangulaci´on de Delaunay en dimensi´on dtiene serios problemas cuando los puntos no est´an en posici´on general, es decir, cuando hay m´as de d+ 1 puntos sobre la misma esfera de dimensi´on d. Incluso cuando los puntos no est´an claramente en posici´on general existen problemas debidos a los errores de redondeo que comete el ordenador al trabajar con n´umeros en coma flotante. Algunas t´ecnicas que resuelven este problema se han propuesto en [20]. 16 Cap´ıtulo 1. Introducci´on y Objetivos 1.1.3. Avance frontal El m´etodo de avance frontal recibe su nombre del hecho de que en cada etapa del proceso existe un frente formado por todos los posibles lados sobre los cuales se puede generar un nuevo elemento. Los generadores de mallas que utilizan el m´etodo de avance frontal han conseguido gran aceptaci´on gracias a su versatilidad y velocidad [11, 41, 49, 51, 60, 65]. La idea b´asica del m´etodo consiste en marchar hacia el espacio que a´un no ha sido mallado, introduciendo un elemento cada vez. La regi´on que separa la parte del espacio que ya ha sido mallada de aquella que a´un no ha sido mallada se denomina frente. La complejidad del algoritmo de avance frontal es de orden O(Nlog(N)), donde Nes el n´umero de elementos. De forma repetitiva el proceso consta de una etapa de generaci´on de un elemento y de otra etapa de actualizaci´on del frente. El primer frente o l´ınea de avance est´a formada por los puntos generados en el contorno exterior del dominio. Para la formaci´on de un elemento son conocidas dos t´ecnicas distintas: 1. Formaci´on de un elemento dada una discretizaci´on del interior del dominio. T´ecnica propuesta por S.H. Lo [50]. Dados dos nodos A y B consecutivos en la l´ınea de avance, es elegido un nodo (C1oC2) de la discretizaci´on interior del dominio o de la l´ınea de avance con el siguiente criterio: AB C C2 1 Figura 1.14: Elecci´on de un nodo cuando el frente avanza 1.1. Algoritmos de generaci´on de mallas 17 a) Para el nodo C1se define: α1=∆ABC1 AB2+BC2 1+CA2 1 (1.2) β1=∆C1BC2 C1B2+BC2 2+C1C2 2 (1.3) γ1=∆AC1C2 AC2 1+C1C2 2+C2A2(1.4) λ1=max(β1, γ1) (1.5) b) Para el nodo C2se define: α2=∆ABC2 AB2+BC2 2+C2A2(1.6) β2=∆C2BC1 C2B2+BC2 1+C1C2 2 =−β1(1.7) γ2=∆AC2C1 AC2 2+C2C2 1+C1A2=−∂1(1.8) λ2=max(β2, γ2) (1.9) c) Se elige C1si: α1λ1≤α2λ2(1.10) En caso contrario ser´ıa seleccionado el nodo C2. 2. Otra t´ecnica b´asica que permite controlar el tama˜no de los elementos y su forma es la siguiente. Dados dos nodos A y B consecutivos en la l´ınea de avance, es elegido un nodo C de la siguiente forma: 1 A B C1 M C6 δδ 1 Figura 1.15: Determinaci´on de puntos equiespaciados 18 Cap´ıtulo 1. Introducci´on y Objetivos a) Elecci´on de un tama˜no de elemento δ. b) Determinar un punto C1a una distancia δ1de A y B con el siguiente criterio: δ1=δsi 0,55AB < δ < 2AB (1.11) δ1= 0,55AB si 0,55AB > δ (1.12) δ1= 2AB si δ > 2AB (1.13) c) Determinaci´on de los puntos C2, C3, C4, C5yC6que estar´an equiespaciados en el segmento que unir´a los puntos C1yM. d) Determinaci´on de los nodos pertenecientes a la l´ınea de avance que est´en situados dentro de un c´ırculo de radio r= 3AB y centro C1. De estos nodos se escogen los que est´en situados en el semiplano delimitado por la recta que une los puntos A y B y que contienen al punto C1. Estos nodos se ordenan seg´un la distancia al punto C1 y se adjunta a una lista de tal forma que el primero ser´a el m´as cercano aplic´andose a todos ellos una penalizaci´on de valor δ1. e) Determinaci´on del punto de conexi´on C. Este ser´a el primer nodo de la lista que satisfaga las dos condiciones siguientes: 1) El interior del tri´angulo ABC no contiene a ninguno de los restantes nodos de la l´ınea de avance. 2) El segmento CM no corta ninguna de las caras existentes en ese momento en la l´ınea de avance. 1.2. El M´etodo de los Elementos Finitos La mayor´ıa de los fen´omenos f´ısicos son formulados por medio de ecuaciones diferenciales en derivadas parciales. Salvo en ciertos casos particulares, la resoluci´on anal´ıtica de estas ecuaciones resulta imposible. El inter´es pr´actico que tiene conocer soluciones aproximadas de estas ecuaciones, y la rapidez en los c´alculos que actualmente proporcionan los ordenadores, han hecho que los m´etodos num´ericos cobren gran importancia. En concreto, el m´etodo de los elementos finitos es uno de los m´as usados para la simulaci´on num´erica de problemas f´ısicos variados, formulados en t´erminos de ecuaciones diferenciales (problemas de c´alculo de temperaturas, desplazamientos, presi´on, velocidades, 1.2. El M´etodo de los Elementos Finitos 19 campos magn´eticos, etc.). La esencia de este m´etodo consiste en el c´alculo de la soluci´on de la ecuaci´on en algunos puntos del dominio de inter´es, denominados nodos. A partir de estos valores se puede calcular la soluci´on en cualquier otro punto mediante el uso de funciones de interpolaci´on. Los c´alculos precedentes requieren, como primer paso, que el dominio est´e discretizado en subdominios de geometr´ıa simple, los elementos. Para una discretizaci´on geom´etrica considerada, podemos utilizar distintos elementos finitos que geom´etricamente sean iguales. En el caso de elementos finitos geom´etricamente distintos, debe prestarse especial atenci´on a verificar la conformidad del mallado, o aceptar la no conformidad en la soluci´on. El conjunto de estos subdominios se denomina malla del dominio. Esta fase del preproceso es muy importante en la medida en que la generaci´on de un mallado para un dominio de geometr´ıa compleja no es una operaci´on simple y adem´as la calidad de la soluci´on num´erica depende del mallado considerado. As´ı, mallados con caracter´ısticas no aceptables repercuten en la no convergencia del m´etodo de los elementos finitos, o bien hacen disminuir la velocidad de convergencia. Tambi´en es cierto, que la calidad de la soluci´on num´erica est´a limitada por la densidad y la calidad de la discretizaci´on. Para problemas con singularidades, dependiendo del tipo de singularidad, existen adem´as elementos finitos especiales que mantienen la misma regularidad caso de que el problema fuera regular. Sin el uso de m´etodos adaptativos, la malla caso de ser irregular, debe atender cualitativamente al comportamiento de estimadores de error locales para reducir el efecto de poluci´on en los elementos cercanos a los que contienen nodos o fronteras en las singularidades. La malla debe acercarse lo m´as posible a algunos requisitos generales, tales como: 1. Adaptarse lo mejor posible a la forma del objeto. En este sentido, las mallas deben adaptarse a las singularidades, haciendo coincidir nodos y contornos de elementos con las singularidades y cargas puntuales, importantes en resoluci´on p-adaptativa, h-adaptativa o h-p-adaptativa. 2. Presentar mayor densidad de elementos en las zonas donde la soluci´on var´ıe m´as r´apidamente. 3. La transici´on de una zona donde los elementos sean grandes a otra donde 20 Cap´ıtulo 1. Introducci´on y Objetivos sean peque˜nos debe ser gradual. Por tanto, las mallas irregulares deben evitar bruscas variaciones de tama˜no entre un elemento y sus adyacentes, siendo de inter´es el dise˜no de criterios flexibles de definici´on de densidades de mallados en el dominio. 4. No debe haber elementos muy degenerados. El sentido preciso de la palabra degenerado depende del tipo concreto de elementos de que se trate. Por ejemplo, si los elementos son tetraedros, ´estos se deben aproximar lo m´as posible al tetraedro regular. As´ı, un tetraedro ser´a degenerado si es muy plano, o muy puntiagudo, o si la diferencia entre la esfera circunscrita y la inscrita es muy grande. Otros requisitos, como la orientaci´on de los elementos, etc., pueden estar motivados por el problema concreto que se pretenda resolver. 1.2.1. Formulaci´on del problema Una clase importante de problemas que aparecen en f´ısica e ingenier´ıa se puede encuadrar en el siguiente marco variacional abstracto: 00Hallar u ∈V tal que :a(u, v) = f(v)∀v∈V00 (1.14) donde Ves un espacio de Hilbert, a:V×V→IR es una forma bilineal continua y el´ıptica, y f:V→IR es una forma lineal y continua. El teorema de Lax-Milgram asegura entonces la existencia y unicidad de la soluci´on. Sin embargo, en la mayor parte de los casos pr´acticos no se puede encontrar la soluci´on exacta del problema anterior; se buscan entonces soluciones aproximadas. La aproximaci´on general de Galerkin consiste en construir subespacios Vh de Vde dimensi´on finita y resolver el siguiente problema aproximado: 00Hallar uh∈Vhtal que :a(uh, vh) = f(vh)∀vh∈Vh00 (1.15) que equivale a la resoluci´on de un sistema algebraico lineal de ecuaciones. Ejemplos t´ıpicos de espacios Vque aparecen en las aplicaciones son H1(Ω), H1 0(Ω), H2(Ω), H2 0(Ω), etc. Los problemas de contorno asociados a ecuaciones 1.2. El M´etodo de los Elementos Finitos 21 en derivadas parciales el´ıpticas de segundo a cuarto orden lineales son ejemplos que se pueden resolver siguiendo el esquema abstracto anterior. El m´etodo de elementos finitos, en su forma m´as sencilla, es un m´etodo espec´ıfico para construir subespacios de dimensi´on finita de V, con el objetivo de hallar soluciones aproximadas de problemas de contorno. El m´etodo comprende la divisi´on del dominio en el cual se encuentra definido nuestro problema en subdominios, denominados elementos finitos, y el uso de formulaciones d´ebiles o variacionales para la obtenci´on de una soluci´on aproximada sobre la colecci´on de elementos finitos. La gran generalidad y variedad de las ideas que subyacen en el m´etodo permite aplicar ´este, con pleno ´exito, en un amplio espectro de problemas en casi todas las ´areas de la ingenier´ıa y de la f´ısica matem´atica. Por otra parte, la aplicabilidad del m´etodo no se limita a los problemas mencionados sino que se extiende a problemas parab´olicos, hiperb´olicos, no lineales, etc. El punto de partida, y a su vez el aspecto m´as caracter´ıstico del m´etodo de elementos finitos, es la subdivisi´on del dominio Ω , en el que est´a planteado el problema a resolver, en subdominios Kmediante, por ejemplo, una triangulaci´on τdel mismo, de modo que se cumplan las siguientes propiedades: 1. ¯ Ω = [ K∈τK 2. ∀K ∈ τ, Kes cerrado y su interior ◦ Kes no vac´ıo 3. Para cada par K1,K2∈τ, ◦ K1T◦ K2=∅ Finalmente, un elemento finito, siguiendo el formalismo introducido por Ciarlet [14], est´a caracterizado por una terna (K,P,Σ) donde i) Kes un conjunto cerrado de IRdde interior no vac´ıo y de frontera ∂K lipschitziana, donde des la dimensi´on de dominio Ω. ii) Pes un espacio de funciones reales definidas sobre K. iii) Σ es un conjunto finito de formas lineales independientes definidas sobre P. 22 Cap´ıtulo 1. Introducci´on y Objetivos 1.2.2. El problema del refinamiento local Adem´as de la generaci´on de la malla (inicial), se puede considerar el problema del refinamiento local. Un refinamiento de mallado excesivamente fuerte, a´un a nivel local, puede producir inestabilidad num´erica en la resoluci´on, debi´endose atender a trabajar en doble precisi´on y a la resoluci´on de sistemas mal condicionados. Como ya se ha dicho, en muchos problemas de ingenier´ıa que simulan situaciones reales es dif´ıcil prever cu´al ser´ıa la malla apropiada. Se pueden realizar varios intentos con el fin de obtener la soluci´on m´as aproximada si se tiene un l´ımite m´aximo en cuanto al tama˜no de los elementos utilizados o sin utilizar un excesivo n´umero de nodos. En un proceso adaptativo hay dos fases: an´alisis del error cometido con la malla anterior y refinamiento de la malla si procede. Los an´alisis del error a posteriori nos permiten calcular para cada elemento una cota del error total cometido -en el caso de estimadores de error-, o bien -en el caso de indicadores de errornos proporcionan un valor por elemento relativo a otros elementos de la malla. En los m´etodos adaptativos, la utilizaci´on de distintos indicadores de error (en problemas el´ıpticos sobre todo) permite plantear estrategias de resoluci´on adaptativas de refinamientos locales para que cada elemento de la malla tenga un error relativo muy similar, de tal forma que el error global sea menor que el deseado. Muchos de los primeros trabajos de an´alisis del error a posteriori fueron debidos a Babuˇska y Rheinboldt [3, 4]. Otros trabajos importantes en este aspecto son los debidos a Kelly, Gago et al. [26]. En esta referencia se ofrecen algunas estrategias sobre el uso de los estimadores de error a la hora de refinar una malla. Zienkiewicz y Zhu se˜nalaban despu´es [105] una estimaci´on del error simplificada que no requer´ıa el c´alculo de saltos de flujo en las caras. Sin embargo el campo del an´alisis de errores esta a´un abierto en muchos problemas en los que no se dispone de estimadores de error, y, adem´as, donde la eficacia de los estimadores propuestos es discutida. Una vez que se ha obtenido una estimaci´on o indicaci´on del error en cada elemento, se pueden seguir varias estrategias para adaptar la malla. Las dos posibilidades b´asicas son el p-refinamiento (o refinamiento en p) y el hrefinamiento (o refinamiento en h). El p-refinamiento consiste en un aumento 1.4. Estrategias de refinamiento 29 Figura 1.19: Diagrama de flujo de un refinamiento adaptativo 30 Cap´ıtulo 1. Introducci´on y Objetivos formar cuatro tri´angulos. Los algoritmos de refinamiento y desrefinamiento del tipo 4-T de Rivara en dimensi´on dos fueron incorporados al c´odigo de Neptuno [21] donde se incluyen algunas estrategias de refinamiento. Por ejemplo se pueden combinar los algoritmos 2-T y 4-T como hacen Ferragut et al. [23]. Sea τuna triangulaci´on y Nel n´umero de elementos de τ. Supongamos definido un cierto indicador de error por elemento, es decir, para cada tri´angulo t∈τse conoce el valor de ese indicador de error: ηt. Con objeto de comparar el indicador de error local con alguna medida del error global, consideremos el n´umero ηopt definido como: η2opt =1 NX t∈τ ηt2(1.16) Entonces se puede seguir la estrategia siguiente: a) si ηt≥γηopt (con γ≈1) el tri´angulo se divide en cuatro elementos. b) en caso contrario, si ηt≥βηopt con β= 0,5, se divide en dos. Con este tipo de subdivisi´on se pretender´ıa alcanzar r´apidamente una malla con distribuci´on uniforme de error. Otra estrategia de refinamiento incorporada al c´odigo de Neptuno es la siguiente: Sea ηmax = m´ax t∈τ(ηt) (1.17) Si ηt≥γηmax, siendo γ∈[0,1] un par´ametro elegido por el usuario antes de ejecutar el programa, refinamos en cuatro tri´angulos. Es decir, esta estrategia nos permite refinar los tri´angulos cuyo indicador de error est´e por encima de un cierto tanto por ciento del indicador m´aximo. L´ogicamente, c´omo se elija un indicador u otro es de capital importancia para el ´exito de una estrategia adaptativa, as´ı como la elecci´on del par´ametro γ. Por ´ultimo, se puede aplicar tambi´en la siguiente estrategia, que no es m´as que una modificaci´on de la anterior. Con la misma definici´on de ηmax que antes, 1.5. Estructuras de datos geom´etricas 31 proporcionamos al programa un γmin y un γmax, ambos entre 0 y 1. Adem´as fijamos el n´umero de nodos aproximado que, como m´aximo, estamos dispuestos a introducir al refinar y sea ´este NNopt. Sea, adem´as, NN el n´umero de nodos en un determinado momento del mallado. Entonces el programa calcula el par´ametro de refinamiento γperteneciente al intervalo [γmin, γmax]⊂[0,1] m´as cercano al cociente NN NNopt . Se consigue de esta forma una mayor libertad en cuanto a los sucesivos refinamientos. Esta variaci´on del par´ametro depende del n´umero de nodos NNopt que es estimado a priori por el usuario. Esto puede ser ´util sobre todo en procesos en los que el n´umero de nodos puede aumentar o disminuir demasiado, es decir, procesos que hayan incorporado el desrefinamiento [66]. En cualquier caso, si nos interesa obtener r´apidos refinamientos, por ejemplo para hallar una malla relativamente fina, cuando se ha partido de una malla inicial grosera, habr´a que elegir valores de γpeque˜nos. Con γ= 0 se obtiene un refinamiento global. Para refinamientos m´as localizados, el par´ametro que habr´a que tomar estar´a m´as cercano a la unidad. En problemas que requieren un alto n´umero de refinamientos y cuyo comportamiento no es conocido a priori, la tercera estrategia, con par´ametro variable, parece la m´as adecuada. Por otra parte, el problema de refinamiento local tambi´en puede extenderse a mallas estructuradas y a otro ´ambito que no sea exclusivamente el de los elementos finitos. Por ejemplo en la figura 1.20 se presenta un ejemplo de h-refinamiento tipo remallado, sobre una malla estructurada quadtree. El ejemplo fue realizado con el software Matlab v.5.0 con la rutina qtdecomp y centra su atenci´on en la detecci´on de bordes de una imagen digital de entrada. Se presentan tres mallas distintas correspondiendo cada una a un valor de 0.83, 0.78 y 0.28 que mide grado de detecci´on deseado (entre 0 y 1). 1.5. Estructuras de datos geom´etricas El dise˜no de estructuras de datos puede considerarse como un ´area de consideraci´on especial en las ciencias de la computaci´on y matem´atica aplicada. La organizaci´on de la informaci´on de cara a la implementaci´on de algoritmos por ordenador exige un esfuerzo necesario de dise˜no que incide claramente en la eficiencia del algoritmo. Teniendo en cuenta que la informaci´on a represen- 32 Cap´ıtulo 1. Introducci´on y Objetivos (a) Imagen digital de entrada (b) Diferentes casos de remallado Figura 1.20: Refinamiento en la detecci´on de bordes 1.5. Estructuras de datos geom´etricas 33 tar (tipo de datos) por las estructuras de datos puede ser de muy diversas caracter´ısticas, como por ejemplo, texto, geom´etrica, gr´afica, etc. tambi´en se exige un esfuerzo de an´alisis asociado al tipo de datos. Es decir, es necesario un estudio de las caracter´ısticas espec´ıficas del tipo de datos, que satisfaga las necesidades del algoritmo que trata con ellas. El problema de organizar y estructurar los datos en forma de estructuas de datos, fue resuelto en la era inicial de la computaci´on de forma simple. La simplicidad era respuesta necesaria ya que los dispositivos de entrada de datos al ordenador eran sencillos. Adem´as, los primeros lenguajes de programaci´on ofrec´ıan posibilidades limitadas en la estructuraci´on y manejo de los datos. Sin embargo, con el avance de la tecnolog´ıa inform´atica y el advenimiento de modernos y sofisticados lenguajes de programaci´on, esta respuesta no es inmediata y ofrece m´ultiples consideraciones y posibilidades. Puede hablarse de que se ha pasado a una sofisticaci´on de las estructuras de datos motivada por la necesidad de dar soluci´on a problemas de creciente complejidad. Los problemas de estructuraci´on de datos que hacen referencia a objetos geom´etricos var´ıan y pueden clasificarse en: 1. Est´aticos. En este caso, todos los objetos geom´etricos en el dominio del problema se proporcionan como parte de la entrada. 2. En-l´ınea (On-line). Donde se permite a˜nadir nuevos objetos geom´etricos pero no permiten ser eliminados hasta que el problema finalice. 3. Din´amicos. Este es el caso m´as general e implica que se permite a˜nadir nuevos objetos as´ı como eliminar otros que ya no hacen falta en la soluci´on del problema. A continuaci´on se enumera una serie de objetivos deseables de las estructuras de datos que tratan con objetos geom´etricos: 1. Capturar la informaci´on estructural (por ejemplo de la malla). 2. Permitir un procesamiento eficiente de las consultas de informaci´on. 3. Permitir actualizaciones de la informaci´on lo m´as eficiente. 34 Cap´ıtulo 1. Introducci´on y Objetivos 4. Optimizar el espacio de almacenamiento. 5. Almacenamiento eficiente en el sentido de minimizar el n´umero de accesos de entrada/salida, cuando el volumen de datos es grande. Concretamente, en el campo de la generaci´on y refinamiento/desrefinamiento de mallas, la necesidad de una buena estructura de datos es crucial debido, entre otras, a las siguientes causas: 1. Entidades de origen espacial (geom´etrico y topol´ogico). 2. Complejidad creciente en la dimensi´on del problema (2D, 3D ...). 3. Volumen de datos considerable, que oscila entre miles y millones de elementos. 4. Creaci´on, eliminaci´on y actualizaci´on de forma din´amica de los elementos. Todo esto ha hecho surgir ´areas espec´ıficas de dise˜no de estructuras de datos. En el campo que nos ocupa se denominan estructuras de datos geom´etricas y estructuras de datos espaciales. Estructuras de datos cl´asicas como vectores, listas, ´arboles o grafos, no son por s´ı mismas muy adecuadas para la representaci´on de objetos geom´etricos, ya que fundamentalmente est´an concebidas para problemas en dimensi´on 1 o bien no son capaces de representar las caracter´ısticas espaciales de los objetos geom´etricos ni del algoritmo que las utiliza. Por ejemplo, en elementos finitos habitualmente es necesario mantener un orden determinado (sentido horario o antihorario) de las aristas dentro de un tri´angulo, o de las caras dentro de un tetraedro. Otras veces, por necesidad de los algoritmos de refinamiento o desrefinamiento, se necesita tener f´acil acceso a los elementos vecinos a uno dado. Este tipo de necesidades se complican al aumentar la dimensi´on del problema, como se ha apuntado anteriormente. De esta forma, es preciso que la combinaci´on adecuada de estructuras de datos o incluso el dise˜no de otras nuevas se adapten a los requerimientos del problema en cuesti´on. Otro aspecto a tener en cuenta es que, por regla general, el dise˜no de estructuras de datos viene estrechamente relacionado con las caracter´ısticas 1.6. Aportaci´on y objetivos de esta tesis 35 del algoritmo que las use. As´ı, un algoritmo que por ejemplo implemente la generaci´on y/o refinamiento de mallas para elementos finitos tipo Delaunay tendr´a unas particularidades en el contexto de la informaci´on geom´etrica a manejar que no coincide con otros algoritmos como el Avance Frontal u Octree. Ello implica que las estructuras de datos de dichos algoritmos no podr´an ser las mismas a cierto nivel de detalle. Por ello se requiere un cuidadoso dise˜no de las estructuras de datos, que adem´as sean lo m´as gen´ericas posibles, pudi´endose as´ı extender su uso a diversos problemas afines. En el cap´ıtulo 3 se hace un estudio y recorrido por las estructuras de datos que en la literatura se usan en el contexto de la generaci´on de mallas y refinamiento/desrefinamiento. 1.6. Aportaci´on y objetivos de esta tesis La aportaci´on fundamental de esta tesis es el desarrollo y estudio de una clase de estructuras de datos geom´etricas basadas en la teor´ıa de grafos. Dichas estructuras de datos surgen de forma natural de un concepto geom´etrico, esqueleto de una malla, que tambi´en es estudiado en esta tesis. De las nuevas estructuras de datos introducidas surgen nuevas versiones de los algoritmos 2DSBR y 3D-SBR, que siguen siendo de complejidad lineal y permiten reducir el costo de espacio y asimismo admiten una m´as clara y completa descripci´on. Objetivos A continuaci´on se presentan los cuatro objetivos fundamentales que se persiguen en esta tesis: 1. Desarrollo de unas estructuras de datos tipo grafo basadas en el esqueleto que de forma eficiente modelen el refinamiento y desrefinamiento por bisecci´on de lado mayor en 2D y 3D. 2. Estudio de las propiedades geom´etricas del grafo basado en el esqueleto y de la partici´on 4T-LE. 36 Cap´ıtulo 1. Introducci´on y Objetivos 3. Implementaci´on de los algoritmos de refinamiento y desrefinamiento basados en el esqueleto en 2D seg´un las nuevas estructuras de datos propuestas. An´alisis de las nuevas versiones de los algoritmos en comparaci´on con los anteriores: 4. Presentar aplicaciones de los algoritmos a: a) Modelos Digitales del Terreno. b) Niveles de Detalle en Visualizaci´on. c) Incorporaci´on de los algoritmos a un c´odigo de elementos finitos comercial y resoluci´on de un problema no lineal. Cap´ıtulo 2 Particiones y algoritmos En algoritmos que tratan con entidades geom´etricas, como los algoritmos de refinamiento y desrefinamiento, es importante estudiar en detalle las transformaciones geom´etricas, en nuestro caso, la partici´on de elementos de una malla. En este cap´ıtulo se hace una revisi´on de las distintas formas de subdividir elementos en una malla y de los algoritmos de refinamiento que surgen de las particiones. Adem´as, se estudian las particiones 4T-LE y 8T-LE as´ı como los algoritmos basados en el esqueleto, por ser de inter´es en esta tesis. Finalmente se procede a estudiar el comportamiento asint´otico de las particiones, proporcionando datos y resultados usados en los restantes cap´ıtulos. 2.1. Definiciones y preliminares Introduciremos las siguientes definiciones y notaciones procedente de geometr´ıa y geometr´ıa computacional que ser´an usadas en adelante. Definici´on 2.1.1 (m-simplex ) Sea V={X0, X1,...,Xm}un conjunto de m+ 1 puntos en IRn(1 ≤m≤n)tal que {−→ X0Xi: 1 ≤i≤m}es un conjunto de vectores linealmente independientes en IRn. Entonces la envoltura convexa cerrada de Vdenotada por S=< V >=< X0, X1,...,Xm>se denominar´a un simplex m-dimensional o un m-simplex en IRn, mientras que los puntos X0,...,Xmse denominar´an v´ertices de S [43]. N´otese que el conjunto Spuede tambi´en definirse de forma expl´ıcita como: S={x∈IRn/x=Pm i=0 λixicon Pm i=0 λi= 1 y 0 ≤λi≤1 para 0 ≤i≤m} 37 38 Cap´ıtulo 2. Particiones y algoritmos A lo largo de esta tesis se usar´an los t´erminos simplex y s´ımplice indistintamente. De esta manera un tetraedro (cerrado) se define por el conjunto (no ordenado) de sus v´ertices < X0, X1, X2, X3>, y una cara por sus tres v´ertices < Xi, Xj, Xk>. Una arista queda definida por un par de v´ertices distintos < Xi, Xj>. V´ertices, aristas, caras y tetraedros son, respectivamente 0−,1−,2−,y 3−s´ımplices. Cualquier i-simplex contenido en un n-simplex S (con i < n) ser´a denominado i-simplex ´o i-cara de S. Como no importa el orden de los v´ertices para definir las caras, resulta que el n´umero de k-caras de un m-s´ımplice es igual a m+1 k+1 . As´ı por ejemplo el n´umero de (m−1)-caras es m+1 m=m+1 1=m+ 1. Definici´on 2.1.2 (Malla conforme de s´ımplices) Sea Ωun conjunto acotado en IRn(n=2 o 3, aunque la definici´on es v´alida en general) con interior no vac´ıo, ◦ Ω6=∅, y de frontera poligonal ∂Ω. Una partici´on de Ωen n−s´ımplices τ={t1,...,tn}tal que (i) Ω = ∪ti; (ii) ◦ ti∩◦ tj=∅si i6=j; (iii) ◦ ti6=∅, se denominar´a una triangulaci´on de Ω. Dos s´ımplices ti, tjde τse dicen adyacentes si ti∩tj6=∅. Si τes tal que s´ımplices adyacentes comparten una cara entera, o una arista entera o un v´ertice, se dice que τes una triangulaci´on conforme o tambi´en que no hay nodos no conformes en τ(en dimensi´on tres se dir´a una teselaci´on conforme o tambi´en una triangulaci´on conforme en dimensi´on tres), ver figura 2.1. Dos triangulaciones (conformes) τyτ∗de un mismo conjunto acotado Ω se llamar´an encajadas oanidadas, y escribiremos τ < τ∗si se cumple la siguiente condici´on: ∀t∈τ, ∃t1, . . . , tp∈τ∗tal que t=t1∪...∪tp. Diremos tambi´en que τes m´as grosera que τ∗, o que τ∗es m´as fina que τ. Definici´on 2.1.3 (Esqueleto) Sea τuna malla de n-s´ımplices. El conjunto skt(τ) = {f:fes una (n−1)-cara de alg´un ti, con ti∈τ}se denominar´a esqueleto o (n−1)-esqueleto de τ, tambi´en denotado por (n−1)-skt(τ) [9]. Por ejemplo, el esqueleto de una triangulaci´on en dimensi´on tres est´a compuesto por las caras triangulares de los tetraedros, y en dimensi´on dos el esqueleto es el conjunto de las aristas de los tri´angulos. Es de destacar por otra 2.2. Algoritmos de refinamiento y particiones en 2D 45 El algoritmo desarrollado por Bank y Sherman [5] (ver figura 2.5 (a)) subdivide un tri´angulo como se ha dicho antes. Se mantiene la regularidad en cuanto a los ´angulos. Por tanto, si realizamos refinamientos globales mediante este algoritmo y representamos por τ0la triangulaci´on inicial y por τnlas obtenidas tras nrefinamientos globales, se tendr´a que α(τ0) = α(τn), con lo que se conserva el ´angulo m´ınimo, puesto que todos los tri´angulos generados son semejantes al original. Sin embargo, si se refina localmente, con objeto de asegurar la continuidad de la soluci´on, Bank introduce la denominada green division (ver figura 2.7) en la que los nodos que se encuentran en el punto medio de alguna cara de alg´un tri´angulo de la malla m´as fina, es decir, los llamados nodos no conformes, se unen al v´ertice opuesto o a otro v´ertice no conforme mediante la inclusi´on de una arista. Esta estrategia tiene el inconveniente de que los tri´angulos as´ı divididos pierden las buenas propiedades de los divididos regularmente respecto de su forma, aunque por otra parte, tiene la ventaja de que la conformidad se logra sin a˜nadir ning´un nodo adicional. Figura 2.7: 00Green division00 para alcanzar la conformidad (a), (b), y (c), seg´un tengamos 1, 2 o 3 nodos no-conformes Definici´on 2.2.2 (Bisecci´on simple y bisecci´on LE) La bisecci´on simple consiste en dividir el tri´angulo en dos sub-tri´angulos uniendo el punto medio de uno de sus lados con su v´ertice opuesto (ver figura 2.5 (b)). Cuando el lado a dividir es el lado mayor del tri´angulo tenemos la bisecci´on por el lado mayor, en lo que sigue, bisecci´on LE (Longest Edge). 46 Cap´ıtulo 2. Particiones y algoritmos Si consideramos la bisecci´on simple, se tiene el problema de que, en general la reiteraci´on de esta forma de refinamiento hace que el ´angulo m´ınimo de la malla tienda r´apidamente a cero. Aparecen tri´angulos degenerados que se comportan mal num´ericamente. En este sentido, la bisecci´on de un tri´angulo se puede mejorar, mediante la bisecci´on por el lado mayor (LE). Se tiene entonces el algoritmo 2-T de Rivara o bisecci´on generalizada. Definici´on 2.2.3 (Partici´on 4T-LE) La partici´on 4-T LE consiste en dividir un tri´angulo en cuatro del siguiente modo (ver figura 2.5 (c)): en primer lugar se divide el tri´angulo inicial mediante bisecci´on por el lado mayor (bisecci´on generalizada) y despu´es mediante segmentos paralelos a los otros dos lados por el punto medio. Rivara [79, 80, 82, 85, 87] ha desarrollado y estudiado un algoritmo eficiente de bisecci´on basado en la partici´on 4T-LE. En la figura 2.6 se presentan las cuatro divisiones que pueden generarse al aplicar un refinamiento basado en la partici´on 4T-LE. Tanto la bisecci´on generalizada, como el algoritmo 4T-LE tienen muy buenas propiedades en cuanto a la regularidad de los elementos a los que da lugar. Un algoritmo equivalente basado en el esqueleto de una triangulaci´on ha sido desarrollado por Plaza y Carey en [71]. La conformidad en los algoritmos 2-T y 4-T se logra del siguiente modo: si alg´un tri´angulo tiene alg´un nodo no conforme en uno de los lados que no sea el mayor se introduce otro nodo en el lado mayor y se une con el v´ertice opuesto; despu´es un segmento paralelo a uno de los lados hace conforme el nodo inicial. Este proceso puede extenderse a otros tri´angulos vecinos. En este caso se dice que el refinamiento se extiende por conformidad. En estos algoritmos, los ´angulos de los tri´angulos para refinamientos de mallas locales est´an uniformemente acotados entre 0 y π. Otro m´etodo de bisecci´on es el presentado por Mitchell [54]. La arista del tri´angulo que va a ser dividida se determina sin ning´un tipo de c´alculo computacional. Solamente se crean cuatro clases de tri´angulos similares y ocho ´angulos distintos. Se satisface tambi´en la importante condici´on de que los ´angulos est´an limitados entre 0 y π. Tambi´en necesita refinamientos adicionales para conseguir la conformidad de la malla. Es de destacar que el algoritmo de Mitchell equivale al 4T de Rivara si la forma en que se elige la arista de bisecci´on (es decir, la arista por la que se realiza la primera bisecci´on), toma siempre 2.3. Algoritmos de refinamiento y particiones en 3D 47 la arista mayor del tri´angulo. Definici´on 2.2.4 (Partici´on Baric´entrica) Para cualquier tri´angulo, la partici´on Baric´entrica consiste en conectar mediante un segmento recto los v´ertices y los puntos medios de los lados con el baricentro del tri´angulo, ver figura 2.5 (d). 2.3. Algoritmos de refinamiento y particiones en 3D En lo que sigue se presentan diversas definiciones que corresponden a particiones de tetraedros en dimensi´on 3. Para una mejor apreciaci´on gr´afica de las particiones se ha optado por representar no s´olo la vista en 3D sino adem´as una vista auxiliar an´aloga a la de B¨ansch [7], que surge al desarrollar las 4 caras del tetraedro en el plano y transformar dichas caras a tri´angulos equil´ateros. Definici´on 2.3.1 (Partici´on Similar) El tetraedro original es dividido en ocho subtetraedros cortando las cuatro esquinas por planos paralelos a la cara opuesta y dividiendo el octaedro interior resultante en cuatro subtetraedros m´as, ver figura 2.8. C A B D D D C A B D Figura 2.8: Partici´on similar en 3D y vista coplanaria de las caras En la divisi´on del octaedro interno en la partici´on similar, es relevante la elecci´on de la arista interna, ya que ello permite optimizar la forma de los cuatro tetraedros que se generan en su interior. Aunque llamaremos a esta partici´on partici´on similar, hay que se˜nalar que s´olo los cuatro tetraedros hijos que est´an en los v´ertices del tetraedro original son semejantes a ´el. 48 Cap´ıtulo 2. Particiones y algoritmos Definici´on 2.3.2 (Partici´on 8T-LE) Para un tetraedro tcon arista mayor ´unica, la partici´on 8T-LE (8 Tetrahedra Longest-Edge) de tse define as´ı: 1. Bisecci´on por la arista mayor de tproduciendo tetraedros t1, t2. 2. Bisecci´on de tipor el punto medio de la ´unica arista de tique es arista mayor de la cara com´un con el tetraedro original t, produciendo tetraedros tij, para i, j = 1,2. 3. Finalmente, bisecci´on de cada tij por el punto medio de la ´unica arista com´un con el tetraedro original t. La figura 2.9 muestra uno de los posibles casos al aplicar la partici´on 8T-LE. M´as adelante, en este mismo cap´ıtulo, se estudiar´an todos los casos posibles. AB C DC AA B D D D Figura 2.9: Caso ilustrativo de la partici´on 8T-LE Definici´on 2.3.3 (Partici´on 3D-Baric´entrica) Para un tetraedro tla partici´on baric´entrica se define as´ı, ver figura 2.10: 1. A˜nadir un nodo nuevo Pen el baricentro de ty en los baricentros de las caras y aristas de t. 2. En cada cara de taplicar la partici´on baric´entrica 2D de la cara. 3. Unir el baricentro Pcon todos los v´ertices de ty con todos los nuevos nodos a˜nadidos anteriormente. Un algoritmo de refinamiento basado en la partici´on similar es el de Bey [10]. Para refinamientos irregulares a la hora de alcanzar la conformidad, Bey considera cuatro patrones que cubren 25 de los 64 posibles, como se ve en la 2.3. Algoritmos de refinamiento y particiones en 3D 49 C A B D D D C A B D Figura 2.10: Partici´on baric´entrica 3D figura 2.11 (n´otese que hay 26= 64 posibilidades, porque s´olo considera el n´umero de nuevos nodos que aparecen en los puntos medios de las aristas, y no su posici´on relativa de acuerdo con la longitud de las mismas). Para el resto de los casos realiza refinamientos globales. Esto implica que en determinados casos puede producirse un efecto domin´o, en el sentido de tener una propagaci´on excesiva del refinamiento a trav´es de la malla. 4 5 6 13 Figura 2.11: Patrones usados por Bey que cubren (4 + 6 + 12 + 3) = 25 posibilidades E. B¨ansch [7] presenta un algoritmo basado en la elecci´on de una arista como arista de refinamiento global en cada tetraedro, pero debe imponer peque˜nas 50 Cap´ıtulo 2. Particiones y algoritmos perturbaciones en las coordenadas de los nodos para evitar ciertas incompatibilidades. ´ El clasifica los tetraedros en dos tipos, denominados rojo ynegro, dependiendo de la posici´on relativa de las aristas de divisi´on y para ello al principio los abre sobre el plano. En la figura 2.12 se muestran los tetraedros tipo rojo ynegro abiertos, la arista de refinamiento global con doble l´ıneas de puntos y las siguientes en el orden de bisecci´on con una simple l´ınea de puntos. 4P P P P 4 4 2 1 P 3 P P 3 4P P P P 4 4 2 1 P (a) Rojo 4P P P P 4 4 2 1 P 3 P P 3 4P P P P 4 4 2 1 P (b) Negro Figura 2.12: Tetraedros abiertos clasificados como rojo y negro Una vez cerrados los tetraedros, ver la figura 2.13, se observan m´as claramente los tipo rojo y negro, aunque los dos tetraedros que dan cuenta de los tetraedros de tipo rojo son equivalentes bajo una rotaci´on de πradianes a trav´es de un eje que pase por los puntos medios de las aristas con nodos P1y P2y la opuesta con nodos P3yP4. 2.3. Algoritmos de refinamiento y particiones en 3D 51 33 4 P P P PP P 2 1 2 1 P4P (a) Tipo rojo 33 4 P P P PP P 2 1 2 1 P4P (b) Tipo negro Figura 2.13: Tetraedros tipo rojo y negro Algoritmos de refinamiento de mallas de tetraedros basados en la bisecci´on son los de Muthukrishnan et al. [59] y Rivara y Levin [83]. Ambos dividen los tetraedros por bisecci´on por el lado mayor. Para alcanzar la conformmidad, tambi´en se aplica bisecci´on por el lado mayor de tetraedros adyacentes. Aunque no se realiza un estudio te´orico-matem´atico sobre la no degeneraci´on de las mallas creadas por medio de este m´etodo, los datos experimentales aportados dan a entender la buena calidad de las mallas obtenidas por medio de este algoritmo, y que incluso se tiene una cierta propiedad de mejoramiento de la malla cuando se parte de una malla inicial muy mala. Recientemente Plaza y Carey [71] han presentado una generalizaci´on del algoritmo 4T de Rivara a tres dimensiones. Esta extensi´on ha sido llamada algoritmo 3D-SBR (3D-Skeleton Based Refinement). Algunas propiedades geom´etricas de este algoritmo han sido descritas por Rivara y Plaza [88]. El algoritmo trabaja primero sobre las caras triangulares de los tetraedros, el 2-esqueleto de la teselaci´on 3D y despu´es subdivide el interior de cada tetraedro de forma congruente con la subdivisi´on del esqueleto. Esta idea podr´ıa ser 52 Cap´ıtulo 2. Particiones y algoritmos aplicada para desarrollar algoritmos semejantes en dimensiones mayores. El algoritmo se puede utilizar para cualquier malla inicial de tetraedros sin ning´un preproceso. Esta es una propiedad importante con relaci´on a otros algoritmos similares. Como en el refinamiento, tambi´en el desrefinamiento, el algoritmo inverso, est´a basado en el 2-esqueleto y utiliza los mismos patrones de elementos usados al refinar. Ambos algoritmos son de complejidad lineal y han sido aplicados a diversos problemas [63, 72]. Liu y Joe [47, 48] presentan un algoritmo similar al de B¨ansch (QLRB: Quality Local Refinement Based algorithm). Ellos clasifican los tetraedros en cuatro tipos y clasifican los tipos de aristas dependiendo del tipo de tetraedro, como se ve en la figura 2.14. La clasificaci´on de los tetraedros en cuatro tipos se realiza de la forma siguiente: 1. DD: Aparecen aristas marcadas 1 y 2 en aristas opuestas. Siendo la arista etiquetada 1 la arista de referencia. 2. DSS1: La arista opuesta a la arista de referencia del tetraedro es etiquetada 3 y las aristas etiquetadas 2 se encuentran en la misma cara. 3. DSS2: Las aristas tipo 2 se encuentran en aristas opuestas. 4. DSS3: Las aristas tipo 2 se encuentran en la misma cara pero una de ellas es opuesta a la arista de referencia. Tambi´en presentan una medida de la calidad de los tetraedros obtenidos. Se prueba tambi´en que el n´umero de clases de semejanza es acotado de forma que las mallas no pueden degenerar. Se han desarrollado otros estudios similares a los de B¨ansch y Liu y Joe. Kossaczk´y [45] desarrolla una aproximaci´on recursiva. Su algoritmo necesita imponer ciertas restricciones y hacer un preproceso a la malla inicial. Su propuesta es equivalente a las de B¨ansch y Liu y Joe en el sentido de que produce las mismas mallas al refinar. Por otro lado Maubach [53] desarrolla un algoritmo para mallas de n-s´ımplices generadas por reflexi´on. Aunque el algoritmo es v´alido para cualquier dimensi´on y el n´umero de casos de semejanza es acotado, 2.4. Las particiones 4T-LE y 8T-LE 53 1 1 1 1 DD 1 2 1 4 3 3 2 1 4 3 3 2 2 3 3 3 2 3 2 DSS1 DSS2 2 1 4 3 3 2 23 3 2 4 3 3 2 33 3 DSS3 Figura 2.14: Clasificaci´on de los tetraedros seg´un Liu y Joe no puede ser aplicado a cualquier malla inicial de tetraedros. Es necesario un cierto refinamiento adicional para evitar ciertas incompatibilidades. Recientemente Mukherjee [2, 57] ha presentado un algoritmo equivalente al de B¨ansch y prueba la equivalencia con el de Maubach. 2.4. Las particiones 4T-LE y 8T-LE Centramos ahora la atenci´on en las particiones 4T-LE y 8T-LE para la descripci´on de los algoritmos de refinamiento en dimensi´on dos y tres. Ambos algoritmos est´an basados en la bisecci´on y tambi´en ambos realizan una clasificaci´on de las aristas que definir´an las sucesivas bisecciones. Teniendo en cuenta la definici´on 2.1.7 y el lema 2.1.2, si se conocen las aristas de bisecci´on, y el orden en que se toman para subdividir el tetraedro, entonces la subdivisi´on est´a perfectamente definida. En este caso el orden est´a basado en la clasificaci´on de las aristas seg´un su longitud. El Algoritmo 4-T desarrollado por Rivara en 1984 [79, 80, 81] como un tipo particular de algoritmo basado en la bisecci´on por la arista mayor, se puede 54 Cap´ıtulo 2. Particiones y algoritmos describir en funci´on de tres patrones de refinamiento que aparecen en la figura 2.15, donde P es el punto medio de la arista mayor. (b) (c)(a) PPP Figura 2.15: Patrones del refinamiento 4-T Un ejemplo de aplicaci´on del algoritmo 4-T para el refinamiento local de un tri´angulo t, aparece en la figura 2.16. En el caso m´as general, el tri´angulo tpuede tener tri´angulos vecinos que comparten las aristas AC y BC y en tal caso, los pasos de propagaci´on por conformidad se tienen que aplicar de forma consistente. El algoritmo, seg´un Rivara [79, 80, 81] que refina un tri´angulo individual t de una triangulaci´on conforme τse expone en el siguiente pseudoc´odigo: El Algoritmo 4-T de Refinamiento {τ, t} /* Se realiza la divisi´on de t*/ Para cada arista ede t, con vecino asociado t?,hace refinamiento-vecino (t?, e) t←t? Mientras tes no conforme, hace encuentra la arista no conforme ede tcon vecino t? Refinamiento-vecino (t?, e) t←t? Fin del mientras Fin del para refinamiento-vecino (t?, e) Si ees la arista mayor de t,entonces realiza la bisecci´on por la arista mayor de t De lo contrario 2.5. Los algoritmos 2D-SBR y 3D-SBR 61 a esta particular implementaci´on del algoritmo por el lado mayor (LE). Un estudio cuidadoso de los posibles patrones de divisi´on con npuntos, producidos por diferentes posiciones relativas de la arista mayor y de las aristas secundarias de t, para n= 1,2,...,6 (lo que incluye la divisi´on global 8T-LE) permite obtener los siguientes resultados [88]: Teorema 2.4.5 El refinamiento en volumen de la malla que se obtiene por la aplicaci´on del 3D-SBR induce el refinamiento de la superficie del 2D-esqueleto de la malla. Adem´as, la superficie de la malla refinada es id´entica a la malla que se obtiene por aplicaci´on del algoritmo de refinamiento 4T a las caras de t. Teorema 2.4.6 Para cualquier malla de tetraedros conforme τ0, despu´es de un n´umero finito de refinamientos locales por medio del 3D-SBR alrededor de alg´un v´ertice P, se obtiene un n´umero finito de tetraedros que comparten el v´ertice P, cuyos ´angulos s´olidos asociados al v´ertice P ya no se vuelven a refinar m´as aunque el proceso de refinamiento contin´ue en torno a P (propiedad fractal).  N´otese que el teorema anterior es la versi´on en 3D de la propiedad an´aloga en 2D (teorema 2.4.3). 2.5. Los algoritmos 2D-SBR y 3D-SBR A partir de la idea de esqueleto de una triangulaci´on es posible formular un algoritmo general de refinamiento de s´ımplices para dimensi´on ncomo sigue: Algoritmo de Refinamiento (τ,n-s´ımplice, partici´on) Para cada n-s´ımplice a refinar, hacer 1) Divisi´on del 1-esqueleto 2) Divisi´on del 2-esqueleto ... n) Reconstrucci´on del n-s´ımplice Fin del para 62 Cap´ıtulo 2. Particiones y algoritmos El esquema anterior acepta como entrada la triangulaci´on τinicial, el ns´ımplice (o varios) a refinar y la partici´on usada para dividir los elementos. Como salida, el algoritmo produce la nueva triangulaci´on τ. El algoritmo divide los esqueletos de la triangulaci´on dada en orden inverso, hasta el (n−1)- esqueleto y finalmente procede con la reconstrucci´on del n-s´ımplice. Este ´ultimo paso precisa adem´as de cierta informaci´on geom´etrica que se deriva de la partici´on, por ejemplo, creaci´on de nodos internos, conexi´on entre caras y aristas etc. En relaci´on con el algoritmo de refinamiento 4T, una versi´on alternativa es el denominando 2D-SBR -algoritmo de refinamiento en 2D basado en el esqueleto-, desarrollado por Plaza y Carey [68, 71]. Dicho algoritmo se ajusta al esquema algor´ıtmico basado en el esqueleto presentado al inicio. B´asicamente realiza el refinamiento en dos pasos consecutivos de la siguiente forma: (1) Identifica y divide las aristas que intervienen a trav´es de todo el proceso de refinamiento; y (2) divide cada uno de los tri´angulos en el proceso de refinamiento como se expone en la figura 2.15 (que depende del n´umero de aristas divididas de cada tri´angulo). El siguiente esquema resume el algoritmo 2D basado en el esqueleto. N´otese que si bien este esquema hace uso de la idea geom´etrica del esqueleto de una triangulaci´on, esto no se refleja en la estructura de datos utilizada. Precisamente ´ese es uno de los objetivos de este trabajo como ya se ha comentado en la secci´on de objetivos 1.6. El Algoritmo de Refinamiento basado en el Esqueleto (τ, t) Encuentra y divide las aristas (τ, t) Divisi´on de los tri´angulos (τ, t) 2.5. Los algoritmos 2D-SBR y 3D-SBR 63 Encuentra y divide las aristas (τ, t) /* Divide las tres aristas de t*/ Para cada arista ede t,hacer Encuentra el vecino t?de tpor la arista e Mientras la arista mayor de e?de t?sea distinta de e,hacer Divide la arista mayor t?de e? t←t? e←e? Encuentra al vecino t?de tpor la arista e Fin del mientras Fin del para Divisi´on de los tri´angulos (τ, t) Para cada tri´angulo t?que tiene al menos una arista dividida, hacer Si el tri´angulo t?tiene una sola arista dividida, entonces Divide t?por bisecci´on de la arista mayor Si tiene dos aristas divididas, entonces Realiza la divisi´on en tres del tri´angulo t? Si no Realiza la divisi´on en cuatro de t? Fin del si Fin del para La figura 2.21 ilustra el uso de los dos pasos principales de este algoritmo para refinar la triangulaci´on (a) de la figura 2.16. En este sentido, remarcamos los siguientes puntos: 1. El procedimiento de bisecci´on por la arista mayor en el 2D-SBR usa impl´ıcitamente el concepto de propagaci´on de la bisecci´on por la arista mayor introducido y discutido en las referencias [85, 88] con el fin de mejorar los algoritmos basados en la bisecci´on por la arista mayor. Este concepto es equivalente al de grafo asociado a una triangulaci´on como se hizo notar en [68], aunque no se utiliza expl´ıcitamente como estructura de datos. 2. El paso de la bisecci´on de la arista asegura la conformidad de la malla 64 Cap´ıtulo 2. Particiones y algoritmos t (a) (b) Figura 2.21: Ejemplo del algoritmo de refinamiento 2D basado en el esqueleto para refinar t refinada. De hecho la divisi´on de un tri´angulo individual se puede realizar en cualquier orden. 3. El algoritmo 2D basado en el esqueleto y el algoritmo 4T, son equivalentes en el sentido de que producen triangulaciones id´enticas. El algoritmo 3D-SBR de Plaza y Carey [68, 69, 71] se puede representar esquem´aticamente diciendo que para cada tetraedro t(que va a ser refinado) en cualquier malla conforme de tetraedros τ, hace: Algoritmo de Refinamiento basado en el esqueleto 3D (τ, t) Encuentra-y-divide-las-aristas Encuentra-y-divide-las-caras Divide-los-tetraedros Se puede destacar aqu´ı que los procedimientos Encuentra-y-divide-las-aristas yEncuentra-y-divide-las-caras act´uan sobre el esqueleto de la malla tridimensional, es decir, sobre la malla de superficie definida por las caras triangulares de los tetraedros. Estos procedimientos realizan las siguientes tareas (aunque no en ese orden exactamente): (1) b´usqueda tri-dimensional de los tetraedros que participan en el refinamiento, (2) bisecci´on de cada una de las aristas y partici´on en 4-Tri´angulos de cada cara del tri´angulo objetivo, y (3) bisecci´on de cada una de las aristas y partici´on parcial de las caras que participan en el refinamiento por conformidad de los tetraedros vecinos. 2.6. Comportamiento asint´otico de las particiones 65 N´otese que, con peque˜nos cambios los dos primeros procedimientos corresponden a la aplicaci´on al esqueleto de la triangulaci´on inicial del algoritmo 2D-SBR. Por otro lado, el procedimiento Divide-los-tetraedros realiza la partici´on en volumen del conjunto de tetraedros cuyas caras han sido refinadas previamente. En el cap´ıtulo 4 se desarrolla m´as extensamente una versi´on del algoritmo 3D-SBR que usa expl´ıcitamente la estructura de datos tipo grafo. 2.6. Comportamiento asint´otico de las particiones Seguimos en esta secci´on las referencias [74, 75]. La pregunta sobre cu´al es el m´aximo n´umero de esferas en tres dimensiones que pueden tocar la esfera unidad, fue causa de una discusi´on entre Isaac Newton y David Gregory. Newton conjetur´o que ese n´umero era 12 mientras que Gregory pensaba que era 13 el n´umero buscado. No fue sino hasta 180 a˜nos m´as tarde cuando Hoppe (1874) demostr´o que Newton estaba en lo cierto [40]. Para un conjunto convexo K, denotemos por N(K) el m´aximo n´umero de copias de Kque pueden tocar a Ksin solapamiento. N(K) se llama n´umero de Newton de K[62]. Se conoce el valor exacto del n´umero de Newton para la esfera unidad en espacios de dimensi´on 2 y 3. Obviamente N(B2) = 6, y como ya queda dicho N(B3) = 12. El problema de encontrar el n´umero de Newton para esferas en espacios de dimensi´on superior es extremadamente dif´ıcil, y tiene su propio inter´es en el campo de la geometr´ıa computacional (ver [62], p´ag. 266 y sg.). En una malla de s´ımplices, se define el grado de un nodo Ncomo el n´umero de nodos adyacentes a N. 2.6.1. Una medida de irregularidad topol´ogica Los conceptos de n´umero de Newton y de grado de un nodo se han utilizado ´ultimamente, para definir una medida de irregularidad topol´ogica para mallas de s´ımplices (tri´angulos en dimensi´on dos, tetraedros en dimensi´on tres) [25, 93]. Es f´acil observar que una triangulaci´on equil´atera tiene en todos sus 66 Cap´ıtulo 2. Particiones y algoritmos nodos interiores grado igual a 6. De ah´ı que se haya definido como medida de regularidad topol´ogica para triangulaciones en dimensi´on dos y tres: t=1 n n X i=0 δi−D(2.1) siendo δiel grado del nodo i,nel n´umero de nodos de la triangulaci´on, y Del n´umero de Newton para la esfera, en el espacio correspondiente [93]. En esta secci´on se estudia el comportamiento asint´otico del grado de cada nodo, en media, cuando el n´umero de refinamientos globales tiende a infinito. Para este estudio se han usado las funciones de generaci´on ligadas a las relaciones de recurrencia que definen el algoritmo. 2.6.2. Comportamiento asint´otico de la media del grado del nodo Las distintas formas de refinar un tri´angulo descritas en la secci´on 2.2 podemos englobarlas dentro de un tipo gen´erico de particiones que quedan definidas por las siguientes ecuaciones de recurrencia: Nn=dNn−1+eEn−1+fTn−1(2.2) En=bEn−1+cTn−1(2.3) Tn=aTn−1(2.4) donde Nn, EnyTnrepresentan el n´umero de nodos, aristas y tri´angulos respectivamente en cada nivel de refinamiento n, con N0,E0yT0dados por la triangulaci´on inicial. Los par´ametros a, b, c, d, e, f configuran un tipo u otro de partici´on, siendo: a: n´umero de tri´angulos que se crean por cada tri´angulo en cada paso. b: n´umero de aristas que se crean por cada arista existente. c: n´umero de aristas que se crean en el interior de cada tri´angulo. d: n´umero de nodos que se crean por cada nodo. e: n´umero de nodos que se crean sobre cada arista existente. f: n´umero de nodos nuevos en el interior de cada tri´angulo. As´ı tenemos que, para el caso de la partici´on similar 4T y la 4T-LE, los par´ametros son: a= 4, b = 2, c = 3, d = 1, e = 1 y f= 0. N´otese que aunque 2.6. Comportamiento asint´otico de las particiones 67 las ecuaciones de recurrencia est´an asociadas a las particiones, no las caracterizan por completo, pues particiones distintas pueden dar lugar a las mismas ecuaciones de recurrencia. Adem´as se puede se˜nalar que en triangulaciones conformes, los par´ametros a, b, c, d, e yfno son independientes entre s´ı. Todas las triangulaciones que aparecen han de cumplir la relaci´on de Euler-Poincar´e. Nos planteamos el siguiente problema: estudiar el l´ımite del promedio del grado de los nodos cuando el n´umero de refinamientos (globales) realizado, n, tiende a infinito para algoritmos basados en la bisecci´on de s´ımplices. Y, m´as en general, el l´ımite de las relaciones de adyacencia en media, cuando el n´umero de refinamientos tiende a infinito. Lema 2.6.1 Sea una triangulaci´on bidimensional con Tntri´angulos, Nnv´ertices yEnaristas, entonces el promedio de tri´angulos por v´ertice viene dado por la f´ormula µn=3Tn Nn (2.5) y el de aristas por nodo: µn=2En Nn .(2.6) Prueba Estudiemos primero el caso del promedio de tri´angulos por nodo. El grado medio de los nodos de la triangulaci´on es 1 Nn Nn P k=1 δk, entonces basta con calcular Nn P k=1 δk, y comprobar que este resultado es 3Tn. En efecto, si realizamos el c´omputo por tri´angulo, el grado de cada v´ertice se incrementa en una unidad, por tanto la suma del grado de los nodos es 3Tn. De igual forma se demuestra para el caso de promedio de aristas por nodo, notando en este caso que cada arista posee dos v´ertices.  Lema 2.6.2 Sea una triangulaci´on tridimensional con Tntetraedros y Nn v´ertices, entonces el promedio de tetraedros por v´ertice viene dado por la f´ormula µn=4Tn Nn (2.7) 68 Cap´ıtulo 2. Particiones y algoritmos Prueba Aplicando el mismo razonamiento de la prueba del lema 2.6.1 se llega a la prueba, notando que en este caso cada tetraedro tiene cuatro v´ertices.  2.6.3. Funciones de generaci´on Una forma de obtener el n´umero de v´ertices y elementos (tri´angulos en dimensi´on dos, tetraedros en dimensi´on tres) para estas particiones, es hallar primero las funciones de generaci´on asociadas a las relaciones de recurrencia que vienen definidas por las propias particiones. La obtenci´on de estos resultados depende de la existencia de estas relaciones de recurrencia. Sean estas funciones N(x), T(x) para los v´ertices y los elementos respectivamente. Una vez halladas las funciones de generaci´on, es f´acil, mediante expansi´on en serie de potencias hallar las expresiones gen´ericas para el n´umero de v´ertices y el n´umero de elementos tras nrefinamientos globales: Nn, y Tn, ver [32, 101], pues ´estos son respectivamente, los coeficientes de los t´erminos de potencia n-´esima para la xen: N(x) = X n≥0 Nnxn, y T(x) = X n≥0 Tnxn, respectivamente. Teorema 2.6.1 Sea τ0una triangulaci´on bidimensional conforme inicial, en la que efectuamos particiones de los tri´angulos basadas en las ecuaciones de recurrencia (2.2), (2.3), (2.4). Entonces el l´ımite del promedio de tri´angulos por nodo tras nrefinamientos globales es igual al n´umero de Newton en dimensi´on 2: l´ım n→∞ 3Tn Nn = 6 Prueba Algunas relaciones entre los par´ametros a, b, c, d, e yfson: 1. d= 1; lo que es obvio ya que en cada paso de refinamiento no se generan nodos nuevos a partir de los existentes. 2. e=b−1; pues si se inserta un nodo en una arista, la arista queda dividida en dos; si se insertan dos nodos, se divide en tres y as´ı sucesivamente. Puesto que el par´ametro eindica el n´umero de nodos que se inserta en cada arista existente y el par´ametro bel n´umero de aristas que se generan por cada arista antigua, se obtiene la relaci´on apuntada. 2.6. Comportamiento asint´otico de las particiones 69 3. f= 1 −a+c; seg´un la f´ormula de Euler: N−E+T= 1 aplicada a las ecuaciones de recurrencia para n= 1, tenemos: T1=aT0 E1=bE0+cT0 N1=dN0+eE0+fT0 Aplicando la f´ormula de Euler y las relaciones anteriores, 1 y 2, se llega a la relaci´on f= 1 −a+c. De las relaciones anteriores, las ecuaciones de recurrencia pueden escribirse como sigue: Nn=Nn−1+ (b−1)En−1+ (1 + a 2−3b 2)Tn−1 En=bEn−1+3(a−b) 2Tn−1 Tn=aTn−1 (2.8) Comenzamos resolviendo la ecuaci´on de recurrencia de los tri´angulos: Tn= aTn−1. Considerando la condici´on inicial, es decir, que inicialmente tenemos T0 tri´angulos, resulta que la funci´on de generaci´on es T(x) = T0 1−ax , o lo que es lo mismo, Tn=T0an. De la expresi´on de T(x) y de la ecuaci´on para las aristas se obtiene: E(x) = bxE(x) + 3(a−b) 2xT(x) + E0(2.9) de donde, sustituyendo T(x) por su valor y despejando E(x), resulta: E(x) = 3 2T0 1−ax + 3 2T0+E0 1−bx (2.10) lo que equivale a que los respectivos valores de Enson: En=3 2T0an+3 2T0+E0bn(2.11) Finalmente, para hallar la funci´on de generaci´on de los nodos, N(x), hemos de utilizar la primera de las ecuaciones de recurrencia, y las funciones ya calculadas E(x) y T(x), resultando: N(x) = xN(x) + (b−1)xE(x) + 1 + a 2−3b 2xT(x) + N0(2.12) 70 Cap´ıtulo 2. Particiones y algoritmos de donde, sustituyendo T(x) y E(x) por sus expresiones y despejando N(x), resulta despu´es de descomponer en fracciones simples: N(x) = T0 2 1−ax + 3 2T0+E0 1−bx +−2T0−E0+N0 1−x(2.13) lo que equivale a que los respectivos valores de Nnson: Nn=T0 2an+3T0 2+E0bn+ 3T0−1 (2.14) Con los resultados que se han obtenido, para calcular el promedio de tri´angulos por nodo, cuando el n´umero de refinamientos tiende a infinito, hacemos: l´ım n→∞ 3Tn Nn = l´ım n→∞ 3anT0 T0 2 = 6 (2.15)  Teorema 2.6.2 Sea τ0una triangulaci´on 3D conforme inicial con N0v´ertices, E0aristas, F0caras, y T0tetraedros. Consideremos la aplicaci´on 8T-LE de forma global. Entonces se produce una secuencia de mallas encajadas τ1, τ2, . . . , τn. En estas condiciones el n´umero medio de tetraedros µnque comparten un mismo v´ertice en la malla τntiende a 24 cuando el n´umero de refinamientos n tiende a infinito. Adem´as, la media del grado de los nodos o media de aristas por nodo tiende a 14. Prueba El refinamiento global de cada elemento de la malla τn−1por la partici´on 8T-LE produce directamente una malla conforme τn(ver teorema 2.4.4). Por las propiedades de la partici´on 8T-LE, los n´umeros de v´ertices, aristas, caras y tetraedros de la malla τn, respectivamente Nn,En,FnyTn, verifican en funci´on de los valores Nn−1,En−1,Fn−1, y Tn−1, las siguientes relaciones de recurrencia: Nn=Nn−1+En−1(2.16) En= 2En−1+ 3Fn−1+Tn−1(2.17) Fn= 4Fn−1+ 8Tn−1(2.18) 3.1. Estructuras de datos en la representaci´on de mallas 77 Lista de aristas: •v1, v2 •v2, v4 •v3, v4 •v4, v5 •v3, v5 •v3, v1 •v1, v6 •v6, v5 •v3, v2 •v6, v3 Una extensi´on a 3D de dichas representaciones es inmediata si consideramos los poliedros que forman la malla en lugar de pol´ıgonos. En este caso es posible definir una lista que considera el conjunto de las caras en funci´on de los v´ertices. La ventaja que tienen estas aproximaciones que definen las entidades de una malla en funci´on de los v´ertices de la misma, es su simplicidad. Si el requerimiento de informaci´on es para la visualizaci´on de la malla, ´estas dos aproximaciones son muy adecuadas, por ejemplo en V RML (Virtual Reality Modelling Languaje) [35], que se est´a convirtiendo en el estandar para la representaci´on y visualizaci´on de objetos 3D as´ı como en OpenInventorT M [98]. Adem´as, la lista de pol´ıgonos ha sido la estructura de datos cl´asica ampliamente utilizada en elementos finitos, y por ejemplo el Partial Differential Equation Toolbox, en Matlab [52] usa este esquema en 2D. Sin embargo, estas estructuras de datos no permiten representaciones completas de mallas ya que no disponen de informaci´on de adyacencia. Aunque tampoco poseen informaci´on de atributos, esto no es un problema ya que es inmediata sin m´as que asignar a cada v´ertice, arista o pol´ıgono (seg´un el caso) dicha informaci´on adicional. Un siguiente paso a la hora de representar una malla 2D es aumentar la informaci´on a almacenar. Consiste en disponer de una primera lista con las coordenadas de v´ertices y dos listas adicionales, una que representa las aristas en funci´on de los v´ertices y la otra que contenga las caras en funci´on de las aristas. Una extensi´on a 3D implica considerar adem´as la lista de los poliedros 78 Cap´ıtulo 3. Estructuras de datos geom´etricas en funci´on de las caras. A continuaci´on se representa el tetraedro de la figura 3.1 (b) atendiendo a esta forma descrita: Lista de aristas + lista de caras + lista de tetraedros. Lista de aristas: •e1:v1, v2 •e2:v1, v3 •e3:v1, v4 •e4:v2, v3 •e5:v2, v4 •e6:v3, v4 Lista de caras: •f1:e1, e2, e4 •f2:e2, e3, e6 •f3:e4, e5, e6 •f4:e1, e3, e5 Lista de tetraedros: •t1:f1, f2, f3, f4 v4   (a) v1 v6 v3 v2 v1 v2 v4 v5 v3 (b) Figura 3.1: Ejemplos de mallas en 2D y 3D 3.1. Estructuras de datos en la representaci´on de mallas 79 Este esquema es el usado por ejemplo en el c´odigo NEPTUNO [21] en 2D y tambi´en en los algoritmos 2D/3D SBR [68, 71] y en el algoritmo de desrefinamiento basado en el esqueleto en 3D [63, 72]. La ventaja de esta representaci´on es que almacena una informaci´on jer´arquica de la malla ya que define las aristas en funci´on de los v´ertices, las caras en funci´on de las aristas y los poliedros en funci´on de las caras. Esta cualidad es muy adecuada para los algoritmos de refinamiento basados en el esqueleto, pues dispone de forma jer´arquica la informaci´on de los conjuntos k-esqueleto para una malla τ. Sin embargo, en l´ıneas generales, esta representaci´on supone almacenamiento excesivo y redundante adem´as de no ser una representaci´on completa por la falta de informaci´on de adyacencia. En las siguientes secciones se estudian estructuras de datos para mallas 2Dy 21 2Dque permiten especificaciones completas de mallas al considerar relaciones de adyacencia entre los elementos de la malla, ver [28, 97]. 3.1.1. Relaciones de adyacencia completas Una forma de representar las relaciones de adyacencia entre las primitivas topol´ogicas del esqueleto es proporcionar enlaces de informaci´on (punteros, desde el punto de vista computacional) a cada una de las primitivas adyacentes a una dada. Seg´un este esquema es posible derivar estructuras de datos basadas en v´ertices, aristas o caras. En la figura 3.2 se presentan esquem´aticamente las tres estructuras de datos derivadas: (a) Basada en v´ertices (b) Basada en aristas (c) Basada en caras. Analizando las tres variantes de la estructura de datos con adyacencia completa, puede notarse que las estructuras de datos basadas en v´ertices y caras figura 3.2 (a) y (c) no poseen tama˜no fijo, ya que el n´umero de punteros depende del n´umero de elementos incidentes al dado. Esta dificultad implica que el coste de almacenamiento es variable y dependiente de la malla en cuesti´on. Asimismo, el coste para hacer una consulta depende de la complejidad local de la malla. En general sucede que las dos estructuras de datos ocupan mucho 80 Cap´ıtulo 3. Estructuras de datos geom´etricas (c) (b) (a) Figura 3.2: Tipos de adyacencia completas espacio y el coste para mantenerlas en todo momento actualizadas es alto. Por contra, la estructura de datos basada en aristas, figura 3.2 (b), no posee estos problemas, ya que a lo m´as, una arista tendr´a dos caras y dos v´ertices adyacentes, salvo cuando pertenece al contorno de una malla, en cuyo caso comparte una. Debido a que la estructura basada en aristas posee esta cualidad la hace eficiente en t´erminos de almacenamiento y se explora en las siguientes secciones, bajo el nombre com´un de estructura de datos basada en aristas, ya que la informaci´on topol´ogica se almacena en relaci´on con las aristas. Para una completa descripci´on de estructuras de datos basadas en aristas ver [97]. 3.1.2. Lista de aristas doblemente conectada (DCEL) Muller y Preparata [58] dise˜naron una estructura de datos denominada Lista de aristas doblemente enlazada, DCEL, (Double Connected Edge List). Esta representaci´on trata cada arista teniendo una direcci´on que impone una orientaci´on a la misma. Avanzaremos alguna definici´on relacionada con grafos en esta secci´on, aunque en el cap´ıtulo 4 se realiza un desarrollo m´as amplio. Un grafo G= (P, L) se dice que est´a contenido en una superficie Scuando se puede dibujar en Stal que dos arcos de Lnunca se cortan. Un grafo es planar si est´a contenido en el plano. Un grafo planar puede siempre incluirse en el plano de forma que todos sus arcos sean segmentos de l´ıneas rectas y que adem´as no se corten entre s´ı. Este tipo de grafos se denomina grafo planar de l´ınea recta, PSLG (Planar 3.1. Estructuras de datos en la representaci´on de mallas 81 Straight Line Graph) o s´ımplemente grafo planar. Los grafos planares juegan un importante papel en geometr´ıa computacional ya que son estructuras que aparecen con cierta frecuencia, por ejemplo en diagramas de Voronoi, triangulaci´on de Delaunay y otras triangulaciones. En lo que sigue se identifican nodos del grafo con nodos de una malla, enlaces del grafo con aristas de la malla, etc. La DCEL para un grafo planar considera cada arco de Lentre (va, vb) conteniendo la siguiente informaci´on: 1. V0, que representa el v´ertice origen (va) 2. Vd, que represente el v´ertice destino (vd) 3. Fl, que representa la cara izquierda seg´un la direcci´on V0aVd 4. Fr, que representa la cara derecha seg´un la direcci´on V0aVd 5. CCW0, que representa el sucesor del arco een el sentido antihorario alrededor de V0, y 6. CCWd, que representa el sucesor del arco een el sentido antihorario alrededor de Vd. e3 f1 f2 ... ... v4v2 3 e f1f0 f2f0 e2v3v1 r FFl v0 vdo d CCW CCW e2 e3v2 f2 1 v f1 e1 e2 e1 v1 v 2 Figura 3.3: Diagrama de la DCEL La figura 3.3 representa un ejemplo de DCEL para el arco e1de la subdivisi´on dada. Seg´un esta estructura de datos, es posible conocer en complejidad lineal las aristas alrededor de una cara dada (en el sentido horario) y las aristas 82 Cap´ıtulo 3. Estructuras de datos geom´etricas alrededor de un v´ertice (en sentido antihorario). Para disponer de esta misma informaci´on pero en otro orden espec´ıfico se deben duplicar los arcos en el grafo planar con orientaciones opuestas a la inicial y de la misma forma construir la DCEL. La complejidad de almacenamiento para esta estructura de datos es O(|V|+|F|+|E|), siendo |V|,|F|,|E|el n´umero de v´ertices, caras y aristas de la malla, respectivamente. 3.1.3. Estructura de datos Winged-Edge Fue desarrollada originalmente por Baumgart [8] a principios de los 70. La estructura de datos Winged-Edge o Arista con Alas se utiliz´o para modelar objetos s´olidos con superficies poligonales en visi´on por ordenador. La potencia de esta estructura de datos hace que sea muy adecuada para el modelado de s´olidos por dos circunstancias: (1) permite modelar el grafo de las adyacencias topol´ogicas entre caras, aristas y v´ertices, y (2) paralelamente a la estructura, existe un conjunto de operadores tipo Euler [38] asociados que posibilita la representaci´on del grafo de las adyacencias topol´ogicas sin utilizar detalles inherentes a la estructura de datos. Dado el grafo planar con v´ertices y aristas, y los pol´ıgonos o caras formando una superficie 2-campo, esta estructura de datos almacena un vector de v´ertices con la informaci´on de la arista adyacente a cada v´ertice. Asimismo, se dispone de un vector de caras que almacena una arista arbitraria adyacente a la cara. Finalmente, la carga de informaci´on se almacena para cada arista e= (va, vb), de ah´ı que puede considerarse esta estructura de datos como basada en aristas, tomando un vector de aristas que almacena los siguientes punteros, ver figura 3.4 : 1. V0, que representa el v´ertice origen (va) 2. Vd, que represente el v´ertice destino (vd) 3. Fl, que representa la cara izquierda seg´un la direcci´on V0aVd 4. Fr, que representa la cara derecha seg´un la direcci´on V0aVd 5. CW0, la arista siguiente a ealrededor de V0en sentido horario 6. CCW0, la arista siguiente a ealrededor de V0en el sentido antihorario 3.1. Estructuras de datos en la representaci´on de mallas 83 7. CWd, la arista siguiente a ealrededor de Vden el sentido horario 8. CCWd, la arista siguiente a ealrededor de Vden el sentido antihorario Arista Adyacente d CWArista Arista CW0 Cara der. Fd Cara izq. Fl Arista CCW0 CCWd Arista Vd Vertice V0 Vertice Arista Adyacente Vertice Arista Cara Figura 3.4: Esquema de almacenamiento Winged-Edge r F e2v3v1f2f0 f0 f1 v4v2 ... ... f1 f2 Fl v0 vdd CCW CW 0CW d e3 e2e4e5 e1 v1 v 2 o 3 e CCW v1 v f2 f1 e2 e4 e5 e3 2 Figura 3.5: Ejemplo de almacenamiento Winged-Edge Para cada arista se almacenan cuatro aristas sucesoras que permiten un acceso eficiente a las aristas alrededor de una cara (sentido antihorario) y de todas las aristas alrededor de un v´ertice (sentido horario). 84 Cap´ıtulo 3. Estructuras de datos geom´etricas Las dos partes sim´etricas en que se puede clasificar cada arista, lo que el autor denomina alas, corresponden a las dos orientaciones posibles para una arista en cuesti´on. Existe sin embargo un caso ineficiente en esta representaci´on que impide el recorrido de forma concisa de las aristas, y es debido a que el puntero a las aristas sucesoras no tiene informaci´on de la orientaci´on concreta que posee la arista. Una extensi´on en este sentido, codifica en un bit de informaci´on la orientaci´on espec´ıfica. La estructura de datos Winged-Edge es similar a la estructura de datos DCEL de la secci´on anterior como se manifiesta en la siguiente proposici´on. Proposici´on 3.1.1 Sea e= (va, vb)una arista cualquiera de una superficie 2-campo, la estructura de datos Winged-Edge es equivalente a la DCEL eliminando los punteros CW0, CWd, la arista siguiente a een torno a V0, Vdrespectivamente en sentido horario. Prueba Se sigue de forma inmediata sin m´as que suprimir los punteros CW0yCWd en la tabla de la figura 3.7 y compararla con la tabla de la figura 3.3.  El coste de almacenamiento para esta estructura es del orden: O(|V|+|F|+ |E|). A continuaci´on se presenta de forma algor´ıtmica los procedimientos que realizan el recorrido de las aristas en torno a una cara en sentido horario, procedimiento WE Arista Cara y el recorrido de las aristas en torno a un v´ertice en sentido horario, WE Arista V ertice. Procedimiento WE Arista Cara (e0) v=e0→v0 e=∅ Mientras e06=ehacer Si v=e→v0entonces v=e→vd e=e→CCWd Si no v=e→v0 3.1. Estructuras de datos en la representaci´on de mallas 85 e=e→CCW0 Fin si Fin mientras Procedimiento WE Arista Vertice (e0, v0) v=e0→v0 e=∅ Mientras e06=ehacer Si v=e→v0entonces e=e→CCW0 v=e→v0 Si no v=e→vd e=e→CCWd Fin si Fin mientras 3.1.4. Estructura de datos Half-Edge La estructura de datos Half-Edge almacena para cada arista en la malla sus dos mitades sim´etricas. En la figura 3.6 se representa una superficie arbitraria resaltando una primera aproximaci´on a la estructura Half-Edge. En esta figura, cabe notar que para cada mitad de arista se almacenan tres punteros, que son: la mitad de arista sim´etrica, la mitad de arista siguiente en el sentido anti-horario y el v´ertice adyacente a la arista. Este primer esquema sin embargo resulta ambiguo ya que no se define expl´ıcitamente ninguna orientaci´on de referencia. Para la estructura de datos Half-Edge, al igual que para las anteriores, es posible asignar una orientaci´on que permite definir de forma ´unica las relaciones con las primitivas adyacentes. De esta forma, asignamos la orientaci´on para cada media arista la que apunta a su v´ertice adyacente superior, denotado simplemente por v. En la figura 3.7 se representa la misma estructura de datos notando los enlaces de informaci´on que se consideran al tener en cuenta la orientaci´on. Puede notarse que para cada mitad de arista se tiene tres tipos de 86 Cap´ıtulo 3. Estructuras de datos geom´etricas Figura 3.6: Esquema de la estructura de datos Half-Edge informaci´on: (1) el v´ertice al final de la media arista, (2) arista que es sim´etrica y opuesta a la actual y (3) arista siguiente en el sentido antihorario. Estas tres informaciones se proporcionan en forma de punteros y en la figura 3.7 se referencian con v,HEnext, y CCW. En total se almacenan 6 punteros para cada arista a representar. Sin embargo, este ´ultimo esquema se puede completar a˜nadiendo m´as informaci´on a cada arista. Concretamente, para cada media arista se almacena su cara adyacente por la derecha y la media arista predecesora en el sentido horario. As´ı se llega a la estructura de datos Half-Edge presentada originalmente por Mantyla [55]. En la figura 3.8 se muestra un ejemplo de este esquema resultante. Se obtiene con este esquema mejorado de la estructura de datos un total de 10 punteros para cada arista a representar. En [44] puede encontrarse una implementaci´on eficiente de dicha estructura de datos para superficies 2-campo. 3.2. Estructuras de datos espaciales 93 3.2. Estructuras de datos espaciales A diferencia de las estructuras de datos goem´etricas estudiadas en la secci´on anterior, en las que el almacenamiento de la informaci´on fundamentalmente es en memoria principal, las estructuras de datos espaciales tienen su mayor inter´es para el almacenamiento y estructuraci´on en memoria secundaria. Cuando la cantidad de datos geom´etricos a manejar es de considerable dimensi´on (varios millones de elementos), el almacenamiento en memoria principal no es eficiente y se recurre al almacenamiento secundario. Los dispositivos de almacenamiento que generalmente considera un ordenador seg´un el modelo Von Neuman, se dividen en almacenamiento primario (memoria principal o RAM) y almacenamiento secundario (disco). Ambos son dispositivos de acceso aleatorio, es decir, una unidad de informaci´on puede ubicarse en cualquier lugar dentro del dominio espacial de la memoria. Los primeros suelen tener capacidades que oscilan en decenas, centenas de megabytes, mientras que los segundos se extienden a varios gigabytes. La memoria principal es m´as cara que los discos. Pero la diferencia fundamental en t´erminos de estructuras de datos viene definida por los dos par´ametros siguientes: 1. Tiempo de acceso. Medido en segundos y representa el tiempo que tarda la CPU de un ordenador en acceder (leer) una unidad de informaci´on del dispositivo. 2. Tama˜no del bloque de transferencia. Mide la capacidad en bits del bloque de datos unidad que se transfiere desde el dispositivo. La siguiente tabla muestra medidas t´ıpicas de estos dos par´ametros en relaci´on a los dispositivos de almacenamiento de memoria principal y disco, tiempo en segundos y tama˜no en bits. En las ´ultimas d´ecadas, el avance tecnol´ogico ha permitido reducir el orden de ambas magnitudes de forma individual. Sin embargo, el ratio entre las dos ha permancecido casi constante. El ratio es importante ya que mide la relaci´on de ambos par´ametros al producirse una transacci´on de informaci´on entre un dispositivo y otro o viceversa. Dicha transacci´on se contabiliza por los accesos a disco y es crucial para las estructuras de datos almacenadas en disco, como 94 Cap´ıtulo 3. Estructuras de datos geom´etricas las estructuras de datos espaciales. Adem´as, se consideran en este tipo de estructuras, dos reglas a tener en cuenta en su dise˜no: 1. Usar una porci´on de almacenamiento principal, que almacena la estructura de los datos ocupados en disco, y 2. Asegurar que la organizaci´on de los datos en disco no degenera cuando se realizan operaciones con los datos de disco, como inserci´on y eliminaci´on. Cabe se˜nalar que aunque las estructuras de datos espaciales que se expondr´an son orientadas a almacenamiento en memoria secundaria, su uso no queda restringido a ´este y es aplicable tambi´en a memoria principal, ver [61, 91]. A continuaci´on se describen las estructuras de datos espaciales m´as importantes. 3.2.1. Acceso Indexado Como idea general, la estructuraci´on por acceso indexado considera el espacio total dividido en dos zonas. En una de ellas, en disco, se almacena los datos en forma no ordenada. En la otra, memoria principal, se construyen unos ´ındices que permiten acceder r´apidamente al bloque en disco donde est´a el dato que se busca. Esta estructura considera b´usqueda por clave ´unica, es decir, a cada unidad de informaci´on almacenada se le asocia una clave que la identifica en el almacenamiento. Concretamente, la zona de ´ındices se almacena como un vector relativo en el que existe un registro por cada posible valor de clave, que se supone num´erica o traducible a num´erica. Cada registro de ´ındice contiene la posici´on en que se encuentra el registro en la zona de datos, la cual tambi´en se maneja como un almcenamiento relativo con los registros ubicados en cualquier sitio libre que hubiese en disco. Esta t´ecnica resulta muy r´apida, pero s´olo es apropiada en el caso de que el rango de valores de la clave no sea muy elevado, ya que de lo contrario el ´ındice Cuadro 3.1: Dispositivos de almacenamiento Par´ametro Memoria principal Disco Ratio Tiempo de acceso 10−7→10−610−2→10−1104→105 Tama˜no de transf. 10 →102104→105102→103 3.2. Estructuras de datos espaciales 95 necesitar´ıa demasiado espacio de almacenamiento. La figura 3.12 muestra un esquema del acceso indexado. Registro c1 Registro c2 Registro cn Registro c4 Registro c3 c1 c2 c3 c4 cn Zona de datos ... ... (memoria principal) (disco) Zona índices Figura 3.12: Esquema de acceso indexado 3.2.2. Acceso Hash La dispersi´on o hashing es una t´ecnica muy utilizada para acceso r´apido a informaci´on almacenada en dispositivos de memoria secundaria. El planteamiento de dos zonas diferenciadas como en el acceso indexado es an´alogo. Los registros en memoria secundaria se estructuran por bloques de celdas, enlazados por listas enlazadas. Existe un vector de celdas en memoria principal con B punteros, cada uno indicando la direcci´on f´ısica para cada bloque en disco. Una funci´on de Hash hhace corresponder cada valor de clave con uno de los enteros [0, B −1]. Si xes una clave, h(x) es el n´umero celda que apunta a una lista enlazada de bloques donde se encuentra la clave pedida x, si tal registro est´a presente en disco, ver figura 3.13. Existen distintas variantes que se puede considerar en el acceso Hash. Primero, si el tama˜no del vector de celdas de la memoria principal es excesivo es posible almacenarlo en memoria secundaria y disponer en memoria principal de otro vector de celdas, que en este caso apunta al vector de celdas de disco que a su vez apunta a otra direci´on de disco en donde se encuentra 96 Cap´ıtulo 3. Estructuras de datos geom´etricas el valor pedido con clave x. As´ı es posible dise˜nar diferentes niveles de acceso Hash. clave x Zona de discoZona memoria principal h(x) Funcion Hash ... ... ... Figura 3.13: Esquema de acceso Hash Por otra parte, es posible hacer que el tama˜no de los registros de disco no sea fijo. As´ı dependiendo de la aplicaci´on en cuesti´on se asigna din´amicamente espacio al registro. Un ejemplo posible se encuentra en el siguiente caso. Imag´ınese que se desea almacenar los tri´angulos que son adyacentes a un v´ertice de una triangulaci´on durante un refinamiento global 4T-LE. Una posible aproximaci´on ser´ıa asignar como tama˜no para 10 tri´angulos adyacentes a un v´ertice. Sin embargo del teorema 2.6.1, p´agina 68 es posible conocer que la media del n´umero de tri´angulos por v´ertice tiende a 6 cuando el n´umero de refinamientos aumenta. Ello induce a pensar que asignar bloques de tama˜no 6 y ampliarlos cuando sea necesario es una idea justificada. Un criterio an´alogo ser´ıa aplicable en 3D, teorema 2.6.2, p´agina 70. En la secci´on 2.6 se encuentra un estudio del comportamiento asint´otico de las particiones del que se derivan resultados aplicables directamente al dise˜no de estructuras de datos tanto geom´etricas como espaciales. Cuando se utiliza este m´etodo el problema principal es la elecci´on de la funci´on Hash. Dado que el n´umero de claves posibles puede ser mayor que el n´umero de ´ındices del vector Hash, habr´a varias claves que se reflejen en un mismo´ındice del vector. Cuando surge este problema se dice que hay colisiones. 3.2. Estructuras de datos espaciales 97 Aunque hay m´etodos de resoluci´on de colisiones, ´estas no son deseables y para evitarlas es necesario que la funci´on que se elija distribuya las claves de forma lo m´as uniforme posible y por tanto se minimice el n´umero de colisiones. Entre las funciones m´as usuales se encuentran: funciones con operador mod, extracci´on de bits de la clave, etc. 3.2.3. Integraci´on arb´orea Un ´arbol es una estructura de datos jer´arquica de una colecci´on de objetos. Se llama grado de un nodo del ´arbol al n´umero de sub´arboles (tambi´en llamados hijos) que tiene. A los nodos con grado cero se les llaman nodos terminales o nodos hojas. El nivel de un nodo es el n´umero de antecesores que tiene desde la ra´ız del ´arbol. La profundidad de un ´arbol se define como el m´aximo de los niveles de los nodos del ´arbol. Se llaman ´arboles generales cuando el n´umero de hijos para un nodo es arbitrario. Un caso particular de ´arboles son los ´arboles binarios, donde cada nodo tiene como m´aximo grado 2. Los ´arboles binarios son ampliamente usados en estructuras de datos y existen diversas variantes tambi´en de uso extendido. A menudo los nodos de un ´arbol binario vienen ordenados de forma que, por ejemplo, los hijos por la izquierda de un nodo dado sean menores que el nodo y los hijos de la derecha mayores. La mayor´ıa de las operaciones como b´usqueda y adici´on de elementos a un ´arbol binario se realizan en O(N), donde Nes el n´umero de nodos en el ´arbol. Diversos autores han aplicado estructuras de datos basadas en ´arboles en alguna parte del proceso de generaci´on y/o refinamiento y desrefinamiento de mallas en 2D y 3D, ver [13, 18, 36, 70]. Su naturaleza jer´arquica permite que se adapten bien a los procesos tambi´en jer´arquicos que surgen en la generaci´on y/o refinamiento y desrefinamiento de mallas. La ampliaci´on multidimensional directa del ´arbol de b´usqueda binario puede orientarse en dos sentidos. El primero de ellos implica la aplicaci´on de la discriminaci´on, en cada nivel, de todas las claves. Entre ellos destacamos el ´arbol cuaternario. La segunda orientaci´on del ´arbol binario a estructura multidimensional hace que en cada nivel se utilice a una de las claves como discriminante. En esta ´ultima categor´ıa destacamos el ´arbol kd. Seguimos acontinuaci´on la referencia [91]. 98 Cap´ıtulo 3. Estructuras de datos geom´etricas ´ Arbol Cuaternario Cada registro tiene asociado un punto del espacio bidimensional que se almacena en un nodo de grado cuatro. La ra´ız del ´arbol divide el espacio en cuatro cuadrantes llamados: NE, NO, SO, SE, por analog´ıa en la orientaci´on espacial geogr´afica. La ra´ız de cada sub´arbol divide a cada uno de los cuadrantes en cuatro subcuadrantes, el proceso se repite recursivamente hasta alcanzar un nodo hoja. Cada nodo consta de cinco campos: un campo [K1, K2] que representa el punto bidimensional, y cuatro m´as que se˜nalan a los cuatro cuadrantes en que se divide el subespacio correspondiente: NE direcci´on noreste, NO direcci´on noroeste, SO direcci´on suroeste y SE direcci´on sureste. En general, para cualquier nodo se enumeran siguiendo el sentido antihorario. En la figura 3.14 se puede apreciar la correspondencia entre un ´arbol cuaternario y los puntos asociados en un espacio 2D. El ´arbol Cuaternario es una estructura que tambi´en puede aplicarse para almacenar la jerarqu´ıa en el refinamiento de mallas triangulares 2D. Concretamente, la partici´on similar 4T y la partici´on 4T-LE, ver definiciones 2.2.1 y 2.2.3 en las p´aginas 43 y 46 respectivamente. En ambos casos es preciso asignar una posici´on de referencia para cada tri´angulo que permita asignar los campos equivalentes NE, NO, SO, SE a los subtri´angulos que se generan en cada partici´on. Para el caso de la partici´on 4T-LE es suficiente tomar como referencia la posici´on que deja el lado mayor a la izquierda cuando se gira en sentido antihorario por las aristas. Un ´arbol cuaternario se dice que est´a optimizado cuando todo nodo cumple la propiedad de que ning´un sub´arbol del nodo puede tener m´as de la mitad de los nodos que el ´arbol cuya ra´ız es nodo. El primer paso para la construcci´on de este ´arbol optimizado es ordenar los puntos por la coordenada xy a continuaci´on por la coordenada y. En el ´arbol resultante la profundidad es del orden de O(log2N). ´ Arbol Kd Un tipo de ´arbol especialmente adecuado para almacenar particiones progresivas de un espacio bidimensional es el ´arbol kd. Los ´arboles kd pueden verse como una implementaci´on eficiente de la generalizaci´on multidimensional de 3.2. Estructuras de datos espaciales 99 B D FC A E G B E D A C F Figura 3.14: ´ Arbol Cuaternario 100 Cap´ıtulo 3. Estructuras de datos geom´etricas los ´arboles binarios y conservan muchas de sus propiedades. En un ´arbol kd, cada registro almacenado en disco se identifica con un nodo del ´arbol. Cada nodo almacena la siguiente informaci´on: 1. Las k-claves K0, ..., Kk−1, que definen un punto del espacio asociado al registro. 2. Dos punteros hizq, hder que son nulos o apuntan hacia otro nodo del ´arbol. 3. Un entero disc, perteneciente al intervalo [0, k−1], que indica el sub´ındice de la clave discriminante, y que en adelante se referenciar´a como discriminante. La condici´on que han de cumplir los nodos de cada uno de los sub´arboles de un nodo se puede expresar como sigue: Si nodo →disc =jentonces Para cualquier nodoien el sub´arbol izquierdo de nodo: nodoi→Kj< nodo →Kj Para cualquier nododen el sub´arbol derecho de nodo: nodod→Kj> nodo →Kj Fin si Algunas reglas para el dicriminante de cada nodo son: 1. Todos los nodos de un nivel dado del ´arbol tienen el mismo discriminante. 2. El nodo ra´ız tiene discriminante 0, sus dos hijos tienen el discriminante 1 y as´ı sucesivamente hasta el nivel k, cuyo discriminante es k-1. El (k+ 1)-´esimo nivel tiene, de nuevo, discriminante 0. De esta forma: Disc →next = (i+ 1) mod k. En la figura 3.15 se muestra un ejemplo de ´arbol kd para representar distintos puntos en una espacio bidimensional de dos claves K0, K1que actuan como coordenadas x, y respectivamente. Las operaciones de b´usqueda tienen en el ´arbol kd un coste O(k·N1−1 k). Para el estudio de otras estructuras arb´oreas avanzadas como el ´arbol bd, ´arbol kdea, ´arbol bm, ´arbol bk ver [61] y [91]. 3.2. Estructuras de datos espaciales 101 K0 K1 G B(10,75) B(40,85) G(10,60) D(25,20) A(50,50) F(70,85) C(80,15) (0,0) (100,0) (100,100)(0,100) A B C D E F Discriminante 0 1 0 1 Figura 3.15: ´ Arbol Kd de dos claves Cap´ıtulo 4 Estructuras de datos tipo grafo 4.1. Definiciones y preliminares La teor´ıa de grafos es particularmente ´util para modelar problemas que manifiestan relaciones entre objetos. En este cap´ıtulo se presenta el grafo del esqueleto de una triangulaci´on en 2D y 3D, que expl´ıcitamente es usado en las estructuras de datos de los algoritmos de refinamiento basados en el esqueleto. Estas estructuras de datos aparecen de forma natural y son consistentes con esta clase de algoritmos. A continuaci´on se recogen algunos conceptos y definiciones de la teor´ıa de grafos usados m´as adelante. Estos pueden encontrarse, por ejemplo, en la referencia [29]. Un grafo geom´etrico es un par (P,L) donde Pes un conjunto no vac´ıo de nodos PyLes un conjunto (posiblemente vac´ıo) de enlaces. El grafo se denota como G(P,L) o s´ımplemente G. Escribimos lk∼(pi,pj)para representar el enlace lkasociado a los nodos (pi, pj). Los nodos representan objetos y los enlaces relaciones entre esos objetos. Los enlaces pueden tener etiquetas opesos y los nodos pueden tener nombres. Un subgrafo G(P0, L0) es un grafo tal que P0est´a contenido en PyL0est´a contenido en L. Un enlace lk∼(pi,pj)entre los nodos piypjse dice que es incidente a los nodos piypj.Los nodos piypj se llaman extremos del enlace (pi,pj). Dos enlaces son adyacentes si tienen un nodo com´un. El grado de un nodo pies el n´umero de enlaces incidentes a pi. Un grafo planar es un grafo que se puede dibujar en el plano conteniendo todos sus 103 110 Cap´ıtulo 4. Estructuras de datos tipo grafo Un aspecto a tener en cuenta es que la adyacencia definida por Tn−1puede restringirse sin p´erdida de generalidad mediante, por ejemplo, la minimizaci´on del n´umero de relaciones que hay que almacenar para un ´unico n-s´ımplice. As´ı que, el siguiente funcional es ´util para este prop´osito, T= m´ın iX j Tn−1 i(j) (4.1) donde el ´ındice irepresenta los diferentes vectores Tn−1para una dimensi´on dada n(n=2,3) y jindexa cada elemento en el vector. En la ecuaci´on (4.1), Tminimiza el n´umero de enlaces en el grafo del esqueleto para un n-s´ımplice. Seg´un esto, las configuraciones de Gn−1para T en dimensiones dos y tres se muestran en la figura 4.4. Grafo del esqueleto 1 T =[1,2,1] 2=[3,1,1,1]T1 T2 2=[2,2,1,1] Vector T Figura 4.4: Grafo del esqueleto y vector T Para nuestro objetivo, el uso del grafo en la modelizaci´on del esquema de refinamiento basado en la bisecci´on por la arista mayor, restringiremos la relaci´on Rentre las k-caras del grafo del (n−1)-esqueleto a aquellas que quedan expl´ıcitamente representadas en la figura 4.4. En la pr´oxima secci´on justificamos esta elecci´on. 4.2. El grafo basado en el esqueleto 111 4.2.1. El grafo del 1-esqueleto en 2D Los algoritmos de refinamiento basados en la bisecci´on por la arista mayor de un tri´angulo en 2D tienen la ventaja de que las diferentes posibilidades de refinar un tri´angulo dependen de determinar la arista mayor. En general, se puede decir que estos algoritmos clasifican las aristas de cada tri´angulo en dos tipos: la arista mayor, de tipo 1, y las otras dos aristas, de tipo 2, y esto es suficiente para nuestros prop´ositos (ver figura 4.5) [71]. Esta idea se puede expresar en t´erminos del grafo del 1-esqueleto al considerar la relaci´on R=“longitud de arista x en el tri´angulo es menor que la longitud de arista y” con el vector T1de la figura 4.4. Como Res una relaci´on de orden entre las aristas, es posible a˜nadir una orientaci´on a los enlaces del grafo que represente este orden. Por tanto, los enlaces en la figura 4.5.(a) indican la dependencia de las aristas al refinar y para asegurar la conformidad. En la figura 4.5.(b) se muestran una sencilla triangulaci´on y el grafo asociado basado en las aristas. Proposici´on 4.2.1 Sea G1el grafo dirigido del 1-esqueleto de una triangulaci´on bidimensional conforme τ. Entonces el grafo no contiene ciclos. Prueba Hay dos situaciones en las pueden ocurrir bucles (ciclos). La primera aparece en mallas que tienen tri´angulos is´osceles o equil´ateros. En el momento de construir el grafo, en los casos en los que hay tri´angulos is´osceles o equil´ateros, se puede dar m´as de una arista de bisecci´on (arista mayor). En estas situaciones en las que se puede elegir m´as de una arista mayor, la primera asignada como arista mayor es la elegida. La segunda posibilidad de bucle es la situaci´on descrita en la figura 4.6. Sea X={x0, x1, x2, x3, ... ,xn}las aristas comunes compartidas por tri´angulos adyacentes en el LEPP de t0(t0es el tri´angulo definido por aristas cyx0). Si aparece un ciclo entonces long(c) <long(x0) <... <long(xn)<long(c). Sin embargo, esto es imposible ya que las aristas de Xson distintas. Si las aristas de Xson iguales, estamos en el primer caso considerado.  Desde el punto de vista computacional la proposici´on 4.2.2 asegura la localidad de la estructura de datos basada en el grafo dirigido del esqueleto. 112 Cap´ıtulo 4. Estructuras de datos tipo grafo AC AB CH I D E GF K J AB CH I D E G J K A B C A 2 C 1 B 2 Arista Tipo B (a) (b) Figura 4.5: Grafo del 1-esqueleto en 2D x3 x0 x1 x2 x3xn c c x0 x1x2 ... xn ? Figura 4.6: Existencia de bucles en el grafo del 1-esqueleto 4.2. El grafo basado en el esqueleto 113 Proposici´on 4.2.2 El coste de actualizaci´on de la estructura de datos basada en el grafo dirigido del 1-esqueleto al aplicar la bisecci´on 4T-LE es como mucho tres nodos en cada paso de refinamiento.  Proposici´on 4.2.3 Sea τuna triangulaci´on conforme bidimensional con N tri´angulos y tun tri´angulo arbitrario de τ. Sea ela arista mayor de t. Entonces el LEPP de tse puede obtener de un recorrido en profundidad del grafo dirigido del 1-esqueleto de τ, comenzando por el nodo edel grafo, con un coste m´aximo del orden de O(N−1). 4.2.2. El grafo del 1y 2-esqueleto en 3D Consideraremos primero el grafo del 1-esqueleto de un tetraedro. Como en 2D, se aplica el mismo criterio en cuanto a la relaci´on entre las aristas. En este caso, se pueden distinguir tres tipos de aristas de acuerdo con su longitud (esta clasificaci´on es semejante a la de B¨ansch [7], Liu y Joe [48] y Mukherjee [57]). La arista mayor del tetraedro la llamamos arista tipo 1. La arista mayor de cada una de las dos caras que no comparten la arista mayor del tetraedro son arista tipo 2, y el resto de aristas son aristas tipo 3 [71]. Observaci´on: Mediante peque˜nas (e hipot´eticas) perturbaciones de las coordenadas de los nodos puede asegurarse que hay exactamente una ´unica arista mayor en cada tetraedro. La arista mayor de un tetraedro es frecuentemente llamada arista de referencia. Para la partici´on 8T-LE, distinguimos tres tipos de tetraedros dependiendo de la posici´on relativa de sus tipos de aristas [71]. Primero consideramos los tetraedros obtenidos al dividir la arista tipo 1, despu´es la o las aristas tipo 2 y as´ı sucesivamente. Cada tipo de tetraedro se puede representar por un grafo orientado del 1-esqueleto como en la figura 4.7. En la figura 4.8 se muestran todas las posibles configuraciones para la partici´on 8T-LE y los grafos orientados del 1-esqueleto asociados. N´otese que en dicha figura se remarca con l´ınea discontinua la arista mayor de cada cara. Observaci´on: Algunas propiedades referentes al grafo dirigido del 1-esqueleto para un tetraedro son: 114 Cap´ıtulo 4. Estructuras de datos tipo grafo 1. El grafo dirigido del 1-esqueleto tiene exactamente 6 nodos y 8 enlaces; 2. El grado del nodo que representa la arista mayor es 4; 3. El grafo dirigido del 1-esqueleto es desconectado y no contiene ciclos. El grafo dirigido del 1-esqueleto de un tetraedro como ha sido presentado anteriormente es una representaci´on basada en las aristas y se puede utilizar para dise˜nar la estructura de datos de cualquier esquema de refinamiento de tetraedros basado en bisecci´on por la arista mayor (ver secci´on 2.2 de la p´agina 43). Esto nos proporciona la informaci´on necesaria para llevar acabo el refinamiento de un tetraedro en t´erminos de sus aristas. Desde el punto de vista computacional, se puede dise˜nar el algoritmo de refinamiento mediante la construcci´on local del grafo para cada tetraedro de forma din´amica. Esta es la idea que usamos aqu´ı por eficiencia en la memoria utilizada. De otra forma, se podr´ıa construir el grafo y almacenar para cada tetraedro en la malla y entonces actualizar el grafo tras cada refinamiento. Esto implica un almacenamiento extra de datos de 8n, siendo nel n´umero de tetraedros en la malla. En la referencia [75] se estudian resultados asint´oticos sobre las relaciones de adyacencia entre los elementos topol´ogicos de las mallas en 3D al aplicar la partici´on 8T-LE. Por ejemplo, cada arista es compartida por 36/7 tetraedros por t´ermino medio. Teniendo esto en cuenta, en el grafo del 1-esqueleto, un nodo podr´ıa tener en media hasta (36/7) ·4 enlaces. Esto implicar´ıa un coste de almacenamiento inaceptable, y por esta raz´on no es eficiente mantener el grafo del 1-esqueleto de forma global a toda la malla, sino localmente para cada tetraedro. 111 3 3223 3 3333 2 33 2 2 (a) (b) (c) Figura 4.7: Patrones del G1para los tres tipos de tetraedros 4.2. El grafo basado en el esqueleto 115 e2 e3 e1 e4 e5 e6 Tipo 1 e2 e3 e1 e4 e5 e6 e6 e5 e4 e2 e3 e1 e4 e5 e6 e2 e3 e1 e4 e5 e6 e6 e5 e4 Tipo 1 e2 e3 e1 e4 e5 e6 e2 e3 e1 e4 e5 e6 e6 e5 e4 Tipo 2 e2 e3 e1 e4 e5 e6 e2 e3 e1 e4 e5 e6 e6 e5 e4 Tipo 3 e2 e3 e1 e4 e5 e6 e2 e3 e1 e4 e5 e6 e6 e5 e4 Tipo 1 e2 e3 e1 e4 e5 e6 e2 e3 e1 e4 e5 e6 e6 e5 e4 Tipo 3  (1) (2) (3) (1) (2) (1) Figura 4.8: Configuraciones posibles de G1para los tres tipos de tetraedros 116 Cap´ıtulo 4. Estructuras de datos tipo grafo En lo que sigue estudiaremos el grafo del 2-esqueleto como una estructura de datos alternativa para modelar una triangulaci´on 3D. Observaci´on: La representaci´on de las figuras 4.9 (a)-(b) es ´unica. En este caso la arista AB es la arista mayor. Esta representaci´on se llama posici´on estandar, ver [7]. Llamamos a las caras f1 yf4 de la figura, caras de referencia. Tenemos dos alternativas para construir el grafo del 2-esqueleto para un ´unico tetraedro, G2 1(τ) y G2 2(τ), ver figura 4.9 (c)-(d). Usando la notaci´on del vector de adyacencia presentado en la definici´on 4.2.2, el grafo G2 1(τ) de la figura 4.9.(c) se escribe como T2 1=[3,1,1,1]. En este caso el nodo con grado 3 es una de las dos caras de referencia (si AB es la arista mayor en el tetraedro, entonces las caras de referencia son f1 yf4). La segunda posibilidad respecto del grafo, G2 2(τ), se muestra en la figura 4.9.(d) y se denota por T2 2=[2,2,1,1]. 4 f f1 f2f3 A D C (a) f4 1 f f2f3 f1 f2 f4 f3 G2 1G2 2 CD D BA D (b) B (c) (d) Figura 4.9: 2-esqueleto para un tetraedro Desde el punto de vista computacional, las configuraciones de los grafos G2 1(τ) y G2 2(τ) son an´alogas para nuestros objetivos. Mediante la construcci´on de dicho grafo, se puede deducir una estructura de datos basada en las caras que es adecuada a las necesidades del algoritmo 3D-SBR como se ve m´as adelante. Proposici´on 4.2.4 Sean G2 i(tj),i, j = 1,2los distintos grafos no dirigidos del 2-esqueleto de dos tetraedros adyacentes t1, t2respectivamente de una triangulaci´on conforme τ. Sea nel nodo com´un ∈G2 i(tj)de las dos componentes de tj, entonces se cumplen las relaciones siguientes para los grados del nodo com´un n. 4.2. El grafo basado en el esqueleto 117 1. Para G2 1(ti), el m´aximo y el m´ınimo del grado de nes 6y4respectivamente. 2. Para G2 2(ti), el m´aximo y el m´ınimo del grado de nes 4y2respectivamente. Prueba Analizaremos las distintas posibilidades de uni´on entre la cara com´un de dos tetraedros adyacentes t1, t2, atendiendo a los grafos G2 1yG2 2. Un nodo que es com´un a las dos componentes de los tetraedros adyacentes t1, t2puede ser una cara de referencia (fref ) para alg´un tetraedro, para los dos o para ninguno. Considerando G2 1, ver figura 4.9.(c), si el nodo en com´un es la cara de referencia para ambos tetraedros, ocurre que el grado de este nodo es 6, y es m´aximo ya que no hay m´as enlaces posibles para este nodo com´un que haga tener un grado mayor, ver figura 4.10 (a) en la que se muestran todas las posibles situaciones en este caso. Si consideramos G2 2, ver figura 4.9.(d), el grado m´aximo es 4, ver figura 4.10.(b). Si el nodo en com´un no es cara de referencia para ning´un tetraedro, el grado m´ınimo de este nodo es 2 para ambos grafos y adem´as es el menor posible para dicho nodo com´un, ver figura 4.10.  Proposici´on 4.2.5 Sea τuna triangulaci´on 3D con ntetraedros y sean, G2 i(tj), i= 1,2yj= 1,...,n, los distintos grafos no dirigidos de los tetraedros de τ. Entonces los grafos son planares y conectados. Prueba Primero veamos que son planares. Recordando la definici´on de grafo planar, un grafo es planar si se puede dibujar en el plano conteniendo todos sus enlaces de forma que los enlaces s´olo se cortan en los nodos del grafo. Una componente de los grafos G2 i(tj) correspondiente a un tetraedro es planar por definici´on, ver figura 4.9. Como los grafos G2 i(tj) est´an basados en las caras de la triangulaci´on, la conexi´on de dos tetraedros es por medio de una cara en com´un y adem´as cada cara es compartida como mucho por dos tetraedros. Las distintas posibilidades de conexi´on entre dos componentes seg´un la proposici´on 4.2.4 son planares y aparecen en la figura 4.10. Entonces, para el grafo completo de una triangulaci´on en la que la conexi´on de componentes es dos a dos tambi´en es planar. Para probar que los grafos son conectados, n´otese que cada componente es conectada ya que una cara es adyacente a dos tetraedros. La 118 Cap´ıtulo 4. Estructuras de datos tipo grafo (fref )fref f f fref f f f fref fref fref fref fref fref (fref ) fref f f fref fref fref fref fref f f f f fff G2 1 G2 2 (b) f f f f f f f f f f f (a) Figura 4.10: Componentes con nodo com´un de G2 1yG2 2 conexi´on sucesiva de componentes conectadas tambi´en origina un grafo que es conectado y esto termina la demostraci´on.  4.3. Implementaci´on del grafo basado en el esqueleto El grafo basado en el esqueleto presentado en la secci´on anterior admite diversas implementaciones, las cuales son comunes a la forma general de representaci´on de grafos, ver [1]. Para representar un grafo dirigido se pueden emplear varias estructuras de datos. La selecci´on apropiada depende de las operaciones que se aplicar´an a los nodos y enlaces del grafo. Una representaci´on com´un para un grafo dirigido G(P, L) es la matriz de adyacencia. La matriz de adyacencia para un grafo G es una matriz Ade dimensi´on n·nde elementos booleanos (uno o cero), donde A[i, j] = 1 si y s´olo si existe un enlace del nodo pial nodo pj. Seg´un esta forma 4.3. Implementaci´on del grafo basado en el esqueleto 119 de representar un grafo, el tiempo de acceso a un elemento es independiente del tama˜no de PyA, y por ello, independiente tambi´en de la malla. La principal desventaja de usar la matriz de adyacencia para representar un grafo dirigido es que requiere un coste de almacenamiento O(n2), incluso si el n´umero de enlaces del grafo es menor que n2. S´olo la operaci´on de leer o examinar la matriz Apuede llevar un tiempo O(n2), lo cual es inaceptable para algoritmos que sean O(n), como por ejemplo los de refinamiento y desrefinamiento basados en el esqueleto en 2D y 3D. En particular, el grafo G1 basado en aristas que es la estructura de datos propuesta para los algoritmos de refinamiento y desrefinamiento basados en la bisecci´on en 2D, tiene la caracter´ıstica de que los enlaces indican direcciones del refinamiento en t´ermino de las aristas, siguiendo el camino LEPP. Un estudio que dejamos para la secci´on 4.6 de la p´agina 133, demuestra que el n´umero de enlaces recorridos al refinar o desrefinar para una triangulaci´on no degenerada, o despu´es de aplicar sucesivos niveles de refinamiento global a una malla es fijo e independiente del n´umero de nodos. Este dato implica que no son muchos los enlaces entre nodos a almacenar para el grafo G1del esqueleto, y que por ello la representaci´on con matriz de adyacencia en s´ı misma no es eficiente. Una mejora notable en la forma de almacenar la matriz de adyacencia es almacenarla como matriz dispersa (matriz sparse). Una matriz sparse es una matriz en la que sus elementos son, en la mayor´ıa, ceros. Dichas matrices son muy comunes en simulaci´on num´erica y surgen cuando la mayor´ıa de las conexiones, en t´ermino de grafos, entre nodos son de caracter local y entre nodos vecinos. Ejemplos de estas matrices son las matrices banda, las matrices diagonales, las matrices triangulares, o simplemente aquellas matrices en la que la mayor´ıa de elementos son ceros. Una matriz sparse se representa por tres vectores de tama˜no q, el n´umero de elementos distintos de ceros, a saber: (1) un vector que almacena los distintos valores distintos de cero, (2) un vector que almacena el ´ındice idel elemento distinto de cero en la matriz original y (2) el ´ındice jdel elemento distinto de cero. Con este nuevo esquema de almacenamiento, una matriz sparse requiere O(q) de almacenamiento, frente a O(n2) con el almacenamiento convencional. Sin embargo, debido a que el acceso en una matriz sparse no es directo, ya 126 Cap´ıtulo 4. Estructuras de datos tipo grafo Algoritmo 3D-SBR (τ, t0) /* Entrada: malla τy tri´angulo a refinar t0 /* Salida: malla τ 1: L= 1 −esqueleto(t0) Para cada arista ei∈Lhacer Subdivision(ei) Fin Para 2: Mientras L6=∅hacer Sea ej∈L Para cada tetraedro ti∈hull(ej)hacer Para cada cara no-conforme fi∈G2(ti)hacer Sea epla arista mayor de fi Subdivision(ep) L=L∪ep Fin para Fin para L=L−ej Fin mientras 3: Para cada cara fi∈G2(τ) a subdividir hacer Subsivision(fi) Fin para 4: Para cada tetraedro ti∈τa subdividir hacer Subsivision(ti) Fin para Fin algoritmo El procedimiento subdivisi´on subdivide los sucesivos esqueletos en orden inverso: pasos 1 y 2 llevan a cabo la subdivisi´on de las aristas, paso 3 realiza la subdivisi´on de las caras y el paso 4 subdivide el interior de los tetraedros. El algoritmo 3D-SBR es de complejidad lineal en el n´umero de nodos [71]. Los puntos donde se usan el grafo del 2-esqueleto G2(τ)se observan en el algoritmo anterior: primero, el bucle interior del paso 2 accede a las caras no conformes de la malla y en el paso 3 en la subdivisi´on del 2-esqueleto. En ambos puntos el grafo del 2-esqueleto proporciona acceso a las caras que toman parte en el refinamiento de forma local y eficiente. 4.5. Algoritmo de desrefinamiento 2D-SBD 127 En el paso 2, el procedimiento hull(e) proporciona el conjunto de s´ımplices S∈τtal que e∈S. Es decir, hull(e) est´a compuesto por todos los s´ımplices que comparten la misma arista e. En [71] se da un algoritmo de complejidad lineal para encontrar la envoltura de una arista (hull). Aunque no se desarrolla implementaci´on del algoritmo 3D-SBR en esta tesis, s´ı puede notarse que la versi´on aqu´ı presentada del algoritmo minimiza el coste de almacenamiento debido a la elecci´on del grafo G2en cualquiera de sus dos opciones, G2 1oG2 2. El grafo almacena la adyacencia de las caras y no es necesario almacenar expl´ıcitamente las aristas de la malla, lo que reduce el coste de almacenamiento de O(|V|+|E|+|F|+|T|) a este otro, O(|V|+|F|+ |T|). En cuanto al coste en tiempo de ejecuci´on, la versi´on aqu´ı presentada no empeora el caracter lineal, ya que en esencia el esquema algor´ıtmico de refinamiento basado en el esqueleto es el mismo. Se deja como trabajo futuro la implementaci´on del algoritmo, en el que presumiblemente la simplicidad y reducci´on de l´ıneas de c´odigo, sean a buen seguro otras caracter´ısticas positivas de la implementaci´on. 4.5. Algoritmo de desrefinamiento 2D-SBD Tras estudiar los algoritmos de refinamiento tanto en el plano como en el espacio en la secci´on anterior, y antes de describir con detalle el algoritmo de desrefinamiento 2D-SBD (Skeleton Based Derefinement) que se ha desarrollado, daremos una serie de definiciones relativas a los mallados encajados no estructurados. Aunque nos limitaremos a tri´angulos, las siguientes definiciones [66, 63] son tambi´en f´acilmente generalizables a otro tipo de elementos y a mayor n´umero de dimensiones. Como ya se dijo en la definici´on 2.1.2 del cap´ıtulo dos, dadas dos triangulaciones de un dominio inicial Ω, τyτ0, diremos que est´an encajadas y que τes m´as fina que τ0, o que τ0es menos fina que τ, y escribiremos τ≥τ0, ´o τ0≤τ, si se verifican las dos condiciones siguientes: 1. todo nodo de τ0pertenece tambi´en a τ. 128 Cap´ıtulo 4. Estructuras de datos tipo grafo 2. toda arista, cara o tetraedro de τ0se puede poner como uni´on (conjuntista) de aristas, caras o tetraedros de τ. Se observa que los respectivos conjuntos de nodos de las triangulaciones verifican la misma relaci´on de contenidos. Es decir, si el conjunto de nodos de la malla τlo representamos por N(τ), se verifica: τ0⊆τ=⇒N(τ0)⊆N(τ) (4.2) La implicaci´on rec´ıproca obviamente es falsa seg´un la definici´on de mallas encajadas. Sea T={τ1≤τ2≤...≤τn}una secuencia de triangulaciones encajadas, yτjuna triangulaci´on cualquiera de T. Se definen: Definici´on 4.5.1 (Nodo propio y heredado) Un nodo Nde τjse dice que es nodo propio de τjsi no pertenece a ning´un nivel de malla anterior, es decir, a ninguna τkcon k < j. En otro caso se dice que el nodo Nes heredado en τj. Evidentemente, por la relaci´on de encajamiento entre los conjuntos de nodos, si un nodo Nes propio de τj, entonces es heredado en τl,∀lcon j < l ≤n Definici´on 4.5.2 (Arista, cara y elementos propios y heredados) Se definen de forma an´aloga a los nodos, es decir una arista, cara o elemento -tetraedro-, se dice propio de un determinado nivel de malla, si no pertenece a ning´un nivel de malla anterior y heredado en ese nivel si fue creado en alg´un nivel anterior. Definici´on 4.5.3 (Aristas, caras y elementos padres e hijos) Se definen de forma natural, por ejemplo, si al refinar, una arista queda dividida en dos, se dice que la primera es la arista padre de estas dos, y ´estas se llaman aristas hijas de aquella. Seg´un el algoritmo de refinamiento que estamos utilizando, una arista solamente podr´a tener a lo sumo dos aristas hijas. Hay que notar tambi´en que no toda arista tiene arista padre. De forma an´aloga se definen las caras padres e hijas. Una cara tiene como mucho, cuatro descendientes directos, es decir, cuatro caras hijas, y, como las aristas, no toda cara triangular tiene cara padre. Por otra parte, una cara solamente puede tener cuatro descendientes directos como m´aximo. Sin embargo, a diferencia de las aristas y de las caras, todo 4.5. Algoritmo de desrefinamiento 2D-SBD 129 elemento, excepto los pertenecientes al primer nivel de malla, tienen elemento padre. Definici´on 4.5.4 (Aristas internas y externas) En dimensi´on dos, al refinar un tri´angulo t∈τjaparecen nuevas aristas, en el nivel de la malla j+ 1. Pues bien, de estas aristas, las que no tienen arista padre se llaman aristas internas, y las que si tienen las llamamos aristas externas. Las primeras se encuentran en el interior del tri´angulo t, las dem´as est´an en su frontera. En dimensi´on tres, definimos las aristas internas y externas a un tetraedro aquellas que pertenecen respectivamente a su interior y a su frontera. Tambi´en denominamos arista interna en volumen como aquella que adem´as de no tener padre, no se encuentran en el interior de una cara triangular, sino que es interna al tetraedro. N´otese, que ´esta es la ´unica nueva arista que no est´a en la frontera del tetraedro que se divide. La figura 4.13 muestra un ejemplo de refinamiento en dimensi´on dos. All´ı, los nodos N1,N2yN3son propios en τjy los nodos A,B,CyDson heredados en τj. Por otro lado, las aristas f1yf2son externas en τj, mientras que f3es interna y c1es heredada. (a) Malla τj−1(b) Malla τj Figura 4.13: Ejemplo de refinamiento Algunas propiedades de las secuencias de triangulaciones anidadas o encajadas en 2D y 3D, son las siguientes: i) Toda arista, cara o elemento propio de cualquier triangulaci´on es heredado en la triangulacion siguiente o tiene sus hijos en ella. Es decir, las familias de aristas, caras y tetraedros relativas a las triangulaciones no est´an encajadas como ocurre con los nodos. 130 Cap´ıtulo 4. Estructuras de datos tipo grafo ii) Por tanto, si un elemento no tiene hijos, pertenece a la malla m´as fina de la cadena de triangulaciones anidadas. iii) De ah´ı que, s´olo aquellos elementos sin sucesores pueden ser eliminados, a la hora de desrefinar, si se quiere salvaguardar la relaci´on de encajamiento de las triangulaciones. Esta ´ultima propiedad es esencial para el algoritmo de desrefinamiento tanto en el plano como en el espacio. Es importante destacar que, como se ver´a m´as adelante, la malla m´as fina va cambiando en el proceso de desrefinamiento y por eso tienen sentido las propiedades anteriores. De esta forma, un elemento perteneciente a un nivel de malla intermedio, ser´a susceptible de ser eliminado al desrefinar, si y s´olo si, previamente, han sido eliminados todos sus sucesores. El criterio de Desrefinamiento En las condiciones de desrefinamiento programadas, se compara la soluci´on num´erica en cada grado de libertad por nodo, con la soluci´on interpolada de la soluci´on en los extremos de su arista entorno (ver definici´on de arista entono en la p´agina 40, definici´on 2.1.5). Hemos utilizado la diferencia relativa y la diferencia absoluta entre estos dos valores. Si esta diferencia es menor que un cierto valor umbral o tolerancia, , que fija el usuario en la entrada de datos del problema, el nodo puede ser eliminado y en caso contrario no. Se puede observar que se comparan dos soluciones aproximadas, no la soluci´on aproximada con la soluci´on exacta, ya que ´esta ´ultima es en general desconocida. Sin embargo, la idea de la condici´on de desrefinamiento es que si la soluci´on interpolada es de por s´ı una buena aproximaci´on de la ´ultima soluci´on num´erica, nos quedamos con la primera. En la figura 4.14, a la izquierda se muestra un nodo propio Ky su arista entorno con nodos extremos en K−1yK+1. Para decidir si se puede eliminar el nodo Ko no, se debe comprobar si para todos los grados de libertad por nodo se cumple la condici´on: ul h(K)−ul i(K)≤(4.3) o la condici´on: 4.5. Algoritmo de desrefinamiento 2D-SBD 131 ul h(K)−ul i(K)≤ul h(K)(4.4) donde ul h(K) es la soluci´on num´erica en el nodo Kcorrespondiente al grado de libertad l, y ul i(K) es la soluci´on interpolada en K, es decir: ul i(K) = ul h(K−1) + ul h(K+1) 2(4.5) Figura 4.14: La condici´on de desrefinamiento La condici´on expresada en la ecuaci´on (4.4) es usada en el ´ambito de los elementos finitos al tratar con soluciones num´ericas. El algoritmo de desrefinamiento 2D-SBD Presentamos el algoritmo de desrefinamiento basado en el esqueleto 2DSBD. En general tenemos una secuencia de mallas encajadas de la forma: T={τ1≤τ2≤... ≤τn}y queremos obtener otra, tras desrefinar T,T0= {τ1≤τ0 2≤...≤τ0 m}, es decir, que se cumple: i) m≤n ii) ∀τ0 k∈T0,∃τi, τj∈Ttales que τi≤τ0 k≤τj Las dos condiciones anteriores se pueden tomar como la definici´on de una relaci´on de encajamiento entre dos secuencias de triangulaciones anidadas definidas sobre el mismo dominio inicial. Es decir que si TyT0son dos secuencias de tales triangulaciones que cumplen las condiciones i) y ii) anteriores, diremos que T0es menos fina que Ty escribiremos T≤T0. Un esquema del algoritmo de desrefinamiento es el que sigue: 132 Cap´ıtulo 4. Estructuras de datos tipo grafo Algoritmo 2D-SBD(T) /* ENTRADA: Secuencia T={τ1< τ2< . . . < τn} /* SALIDA: Secuencia Tm={τ1< τ0 2< . . . < τ0 m} 1: Para j=nhasta 2hacer Para cada nodo propio elegible N∈τj,hacer /* Se calcula el conjunto L∈1-esqueleto con las /* aristas que han de refinarse en el paso 2 Se examina la condici´on de desrefinamiento Se marcan los nodos a eliminar 1.1 Se comprueba la conformidad local 1.2 Si un nodo no se elimina, tampoco los de la arista entorno Fin para Fin para 2: Para j= 2 hasta nhacer /* Se redefine la malla τj τj= 2D−SBR(τj−1, L) Fin para Fin algoritmo El algoritmo de desrefinamiento desarrollado en esta tesis, aunque es an´alogo al de [66] en el sentido de que ante una entrada espec´ıfica produce la misma salida al final, tiene diversas caracter´ısticas que lo diferencia de la versi´on anterior. A continuaci´on presentamos la siguiente comparativa entre los dos algoritmos: 1. La estructura de datos para almacenar la malla pasa de ser de (lista de nodos + lista de aristas en funci´on de los nodos + lista de tri´angulos en funci´on de las aristas) a esta otra, considerada como ´optima, (lista de nodos + lista de tri´angulos en funci´on de nodos). 2. La estructura de datos usada para implementar la partici´on 4T-LE en el algoritmo usa la estructura de datos del grafo G1. 3. Se reducen las l´ıneas de c´odigo. La nueva versi´on implementada en Matlab reduce considerablemente el n´umero de l´ıneas de c´odigo de la versi´on de Plaza, la cual es c´odigo Fortran. 4. La simplicidad de la nueva versi´on queda manifiesta, aparte de por el reducido n´umero de l´ıneas de c´odigo, por la flexibilidad de introducirlo 4.6. Propiedades geom´etricas de la partici´on 4T-LE 133 en diversas aplicaciones. Esta caracter´ıstica es en parte debido a la modularidad que ofrece el lenguaje Matlab. Por ejemplo, se incorpora en el lenguaje de realidad virtual V RML para generar niveles de detalle, ver aplicaciones. 5. El nuevo algoritmo cambia la forma de proceder respecto al anterior. Se puede decir que la versi´on aqu´ı presentada act´ua de forma forward, ya que parte del nivel de malla m´as grosera hasta obtener la malla m´as fina. Esto puede verse como un proceso de reconstrucci´on de la malla en el que el refinamiento es selectivo y de acuerdo con el concepto de desrefinamiento de la malla. Por contra, el algoritmo de [66] es tipo backward, ya que procede a la inversa, desde el nivel de mallas m´as fino a la m´as grosera. 6. La portabilidad del nuevo algoritmo se garantiza debido a las posibilidades de traducci´on de c´odigo Matlab a por ejemplo, lenguaje C o Java. 7. La complejidad lineal del algoritmo en t´erminos de tiempo permanece, sin embargo, la complejidad en cuanto al almacenamiento se reduce y se mejora, ya que no se almacena la lista de aristas de la malla. Para ilustrar el uso del algoritmo de desrefinamiento 2D-SBD, presentamos un ejemplo que consiste en el desrefinamiento de 4 niveles de mallas encajadas. El criterio de desrefinamiento es por diferencia absoluta en los nodos con tolerancia de 0,1875. El ejemplo de desrefinamiento aproxima la singularidad de una funci´on matem´atica en el punto (−1,0). En las las figuras 4.15 y 4.16 se ilustran los cuatro niveles de mallas antes de desrefinar y despu´es de desrefinar (el s´ımbolo asterisco indica que el nodo es eliminado en ese nivel de desrefinamiento). Para completar en dimensi´on 3 con un algoritmo de desrefinamiento y aplicaciones, emplazamos a las referencias [63, 72, 73]. 4.6. Propiedades geom´etricas de la partici´on 4T-LE En esta secci´on se estudian algunas propiedades geom´etricas interesantes de la partici´on 4T-LE para mallas bidimensionales triangulares. Los algoritmos 134 Cap´ıtulo 4. Estructuras de datos tipo grafo −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 (a) Malla 1, p=9, t=8 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 (b) Malla desrefinada, p=9, t=8 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 * * * * (c) Malla 2, p=25, t=32 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 (d) Malla desrefinada, p=21, t=28 Figura 4.15: Mallas 1 y 2 del test del algoritmo 2D-SBD 4.6. Propiedades geom´etricas de la partici´on 4T-LE 135 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 * * * * ** * ** * * * * * ** * * * * * * * * * * * * * * * (a) Malla 3, p=81, t=128 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 (b) Malla desrefinada, p=46, t=76 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 * * * * ** * * * * * * ** * * ** * * * * * * ** * * * * * * * * * * *** ** * * ** * * * ** ** * ** * * ** ** *** * * * *** * * * ** ** * ** * * * ** ** ** * ** ** * * ** * * * *** ** * * * * * * ** * * ** * * * * *** * *** * * * * ** * * (c) Malla 4, p=289, t=512 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 (d) Malla desrefinada, p=131, t=242 Figura 4.16: Mallas 3 y 4 del test del algoritmo 2D-SBD A.3. Communications in Numerical Methods in Engineering 243 ðcñVòmó ôQõ ònö÷ øQøVòùòöù4ñúmû ù4úmñV÷ö)üýmñQöþc÷2ÿ÷2ùmý  Qõ ònö÷ øñV÷2ü  4÷  P÷  Vùhòmÿ*ðnýcñ  Zù4ô  ö   !#"%$'&#(*)+ ,-"./"0&12!#),3.4 '56798;:-<>= <>= ?@9AB9A-6DC2EF?:G:>6IH9@= @J#<>596<>:>= 8;@JBK 8L<>= ?@M= @N?B9:G8K J?:>= <>5OPA8:>6(<>56Q' 4877:>?R8ST5EF?: <>56U;VSW8A-628@9CX<>596NY' Z8;779:>?8S>5XEF?:#M9[Z= <>5/<>56NA-\6IK 6I<>?@]J:T875_^<>596M8;K J?:>= <>59OPA8CON= < 8MOU?:>6`@98;<>B:T8;K_8@9CaSI?@9A-= A3<>6W@<C6DA>Sb:>= 7<>= ?@dcFA-6W62+= JB:>6Reb '56`3$')Z8K J?:>= <>5Of8A(7:>6DA-6I@R<>6DCNghUijK 8;kD88@9CXl8:>6Ih]m nWoLprq'?:>\AG= @P<3q?MA-6DsB96I@R<>= 8;KrA3<T8J6DA 8A E?K K ?;qAWtc=e = C6I@R<>= EFhN8;@rCMg= A-6DSb<G<>566DCJ6DAA-76WSI= H96DCNgRh<>56?Lu6W:T8;K K:>6IH9@6WOU6W@<j= @9C= SW8L<>?:TAW^cFA3<>6W79A n8;@9CN#= @U+= JB:>6eb^;8@9Cac= = eGA-BgC= u= C6<>56= @rC= u= CB98;K<>:>= 8@JK 6DA <>?C6IH9@6'<>56@6Wqv<>:>= 8;@JBK 8L<>= ?@_^ cFA3<>6I7]Rebrl(?@9SI6I:>@9= @JNSI?OU7B<T8;<>= ?@98;K<>= OU6^= <#A-5?BK Cwgr6`@9?;<>6DCw<>598L<<>5623$')08;K J?:>= <T5Ox= A?E K = @6D8;:Sb?OU7K 6by= <-hwq= <>5]:>6DA-7r6DST<<>?U<>56@RB9OMgr6W:?;EG@?C6DAW9"BOU6W:>= SI8K_6Iy7r6W:>= OU6W@<TASW8;:>:>= 6DCw?B<= @ <>5= A:>6W7r?:-<:>6Wu6D8;K9<T598L<<>5= ASI?A3<8;797:>?8ST56DA8A-hROU7<>?<>= SW8K K hP8UA-OU8K KSb?@9A3<T8@<8A'A>B9ST5X8NOU6DA-5X= A :>6bHr@6DC Algorithm 2D-SBR(τ, t0) /* Input: τ mesh, t 0 triangle /* Output: τ mesh 1. L=1-skeleton(t 0) For each edge ei ∈ L do Subdivision(ei) End 2.For each edge ei ∈ L do Si=LEPP(ei,G1 (τ)) For each triangle ti ∈ Si do Let ei be the longest edge of ti Subdivision(ei) End End 3.For each triangle tj ∈ τ to be subdivided do Subdivision(tj) End Algorithm 3D-SBR(τ,t0) /* Input: τ mesh, t 0 tetrahedron /* Output: τ mesh 1. L=1-skeleton(t 0) For each edge ej ∈ L do Subdivision(ej) End 2.While L = ∅ do Let ek be an element ∈ L For each tetrahedron t i ∈ hull(ek) do For each non-conforming face fi ∈ G2(t i) do Let ep be the longest edge of fi Subdivision(ep) L=L∪ep End End L=L-ek End 3.For each face fi ∈ G2(τ) to be subdivided do Subdivision(fi) End 4.For each tetrahedron ti ∈ τ to be subdivided do Subdivision(ti) End z{ |W}~3RG#b'INI' ( |WW~3{ 3RT ,@U<>56SW8A-6?;Er<T563$)Z8;K J?:>= <>59O]^D<>56EF?B9:jOP8;3?:(A3<>6W79A8;:>6t9c3nLeGBgC= uR= A-= ?@P?;E6DCJ6DA<>5r8L< E?:>Ox<>56UnIA-\6IK 6I<>?@X?;E <>596`= @97B<<T6b<>:T856DC:>?@_ce'RBgC= u= A-= ?@X?E 6DCJ6DACB96#<>?N<>56`7:>?798J8;<>= ?@ ?;E<>596aSI?@E?:>OU= <-h<>?<>596w@6W= J5gr?:M<>6I<>:T856DC:T8'cFe2#8Du= @9Jg= A-6DST<>6DC<>596X6DCJ6DA2= OU7K = 6DCv= @v<>56 :>6bHr@6WON6W@R<D^A-BgC= u= C6P<>596XSI?:>:>6WA-7r?@9C= @J]EF8SI6WAM8SWSb?:TC= @Ja<T?]<>56aQ' 798;:-<>= <>= ?@EF?:UMcQRe i6W:-EF?:>O*<T56`8;797:>?79:>= 8;<>6= @<>6W:>= ?:A-BgC= uR= A-= ?@X?;E_<>56#<>6b<T:T8;56DC:>?@U<T8;\R= @JN798:-<'= @w<>56#:>6IH9@96IOU6W@<L^ E?K K ?Lq= @9Jv<>56Y 798;:-<>= <>= ?@0A3<>:T8;<>6WJh,</A-5?BK CZgr67r?= @<>6DC?B<a56W:>6<>598L<a<>567:>?SI6DCB:>6 TrTL¡¢W¡ b¡F£;¤ = @Z+= JB:>6/A-BgC= uR= C6DAU<>56A-B9SWSb6DA>A-= u6]A-\6WK 6I<>?@9AN= @:>6Wu6W:TA-6X?:TC6W:Dt<>6W79A]n]8;@rCZ 7r6I:-EF?:>O¥<>56MA-BgC= uR= A-= ?@?;EG<>5626DCJ6DAW^9A3<T6I7P7r6I:-EF?:>OPA<>56MA-BgC= u= A-= ?@?;EG<>56E8SI6DA8@9C]A3<>6W7Q û r¦b§I¨D©ª «b¬Iw® ¯V°b±b±b±² ¦b¬D³#´Xª µ ¶¨#· Jö I¦b³L¸-¹ ÿ ºR» ¼9½b¾¾'¿WÀÁÂ¿D¾Ã-Ä-Á Å2Ã-Æ ÇLÁ_ÈÀWÉIÀWÉ °b±T±b±DÊË;ËÌ ÍÎWÏ ÐÄFÃFÑDÒbÄFÃÓ#¿DÔÕ ÀWÉÖ3×IØ(ÙbÚWÛ ÜDÝ Ö3Þ ß 244 Ap´endice A. Publicaciones à áDâLãâäWå]æ çGè_éêë;ì â í_â;î çGè_é9ï4çGðñ/ç âDãò çjê;ç órôIõ-öF÷õ>øPùú>ûôPù-üýþÿ  Rÿ ù-ÿ ÷  ÷öjúTôbú>õ  ûôDþõ  'ûô   ÷õ>ÿ ú>ûø ÿ ù#÷;ö  ÿ  9ô  õ  I÷øUó  ô  ÿ ú  ]ÿ  /úTûô  RüøMýrôWõ ÷;ö  ÷þôDù  ù2ù3ú  ;ú>ôDþ4ÿ   'ûôPór÷ÿ  RúTù ! ûôWõ>ôUú>ûô #"%$ ù '& ô  ôIú>÷ ( õ  ;ó9û *)!+-,/.0 ÿ ùü9ù-ôDþ * õ>ô  ô  ;õ 1  ù-ôIô  vÿ  vú>ûôaóõ>ô  Rÿ ÷ürù 23  ÷õ>ÿ ú>û9ø #4 ÿ õTù3ú % ú>ûôXÿ  ôWõ  ÷÷óv÷öù3ú>ôWó 5"( IôWù>ù-ôDùú>ûô # ÷ 67$ I÷ 6 ö÷õ>øUÿ  ö 86 bôDùÿ  ]ú>ûô2øUôDù-û 93 rþaÿ  /ù3ú>ôWó (: Nú>ûôMù-üýþÿ  Rÿ ù-ÿ ÷  /÷;öGú>ûô ;"%$ ù '& ô  ôIú>÷  aÿ ùórôWõ-ö÷õ>øUôWþ <= ]ýr÷ú>û/ù3ú>ôWórù  ú>ûô ;"$ 3ù '& ô  ôbú>÷ 6# õ  ;ó9ûXó9õ>÷ % ÿ þôDù >6 IôWù>ùú>÷Uú>ûô2ö 86 IôWùú 3& ÿ  wó  õ-úÿ  aúTûôMõ>ô ? ôWøUô  Rúÿ ( ÷ 73 9þ ô A@B Iÿ ô  Rú C D%6E= dù3ú>ôWó F"7 ú>ûôwóõ>÷  IôWþüõ>ô BGHEI8I1,8J%0 óõ>÷ % ÿ þôWù`ú>ûôwù-ôbúM÷;öù-ÿ øUó  ÿ  bôDù 2KML5. ù-ü  Tû4ú>û  Lú JCLNK 'û  ;úÿ ù GHEI8I1,8J30 jÿ ù  I÷øUóõ>ÿ ù-ôDþX÷;ö O3  ú>ûô2ù-ÿ øUó  ÿ  IôDùú>û  Lúù-û  ;õ>ô#ú>û9ô2ù 1 ;øUô`ôDþ 7 ô`ô 6 P RQTSUVQW=X(QVYZ>[N\ZY]^_Q[7`!QaQVZ\ZY$bDcZ[QVd]c_aeC\ZW='^_Xf,g"dMhDcZ[7QD0 i ô C I÷  9ù-ÿ þôWõ > Mó  ;õ-ú>ÿ  Iü   õþ÷ø j ;ÿ N b÷õ>õ>ôDù-ór÷  rþÿ  2ú>÷Nú>ûô 2e #õ NhD3 õ>ÿ C= -ù '  9þ *,g[ Ró  ÿ E0D3 rþ B ÷ 73 õ>ô A? ôWøNô  Rúü9ù-ÿ  P R ÷  ôDù3úôDþ 7 ô`ýÿ ù-ô  bú>ÿ ÷  /÷  aú>ûõ>ôWô C I÷  Rú>ÿ  ü÷ü9ùù-üýþ÷ø j ÿ  9ù  i ôö÷  bürù÷  aú>ûô ö÷ 6  ÷ % ÿ  ú  '÷ór÷ÿ  úbù k_,0c þôWøU÷ 6 9ù3ú>õ  Lú>ÿ ÷ 6 ÷;ö#ú>û9ô/ù '& ô  ôIú>÷  ý  ù-ôDþ 53  ÷õ>ÿ ú>ûø l ùPþôDù 1 Iõ>ÿ ýrôDþ ÿ  ú>ûôXó9õ>ô  IôDþÿ  ù-ô  bú>ÿ ÷  rù2÷öú>û9ÿ ùNó  órôWõ m,g"60;[ R÷øUôaù3ú  Lú>ÿ ù3ú>ÿ 3 øNô  ù-üõ>ôDù2÷;öúTûô #noCp!p ù-üý  õ  óû þôIú>ôIõ>øUÿ  ôDþ B Lúô  TûPõ>ô ? ôWøUô  Rúù3ú>ôWó 6= Pú>ûô >? 9õTù3ú'ù3ú>ü9þ ; 'ô _ I÷ 6 9ù-ÿ þôIõjú>ûô !e #õ 3#hD3 õ>ÿ  þ÷ø j ;ÿ q[ ÿ ú>ûNøUôDù-ûNõ>ô ? ôWøUô  új÷  Nþÿ ù 8r 3÷ÿ  úù-üýõ>ô  ÿ ÷ 6 9ù OKs6K +  9þ Kut ÿ ú>û BK*vRKs-wZK + w_Kut6 DöF÷õjÿ  ôWõ>øU÷ù3ú õ>ô  ÿ ÷ 6BKut3 ÿ  Rú>ôIõTøNôDþÿ  LúTôõ>ô  ÿ ÷ BK +  9þU÷üú>ôIõ>øU÷Rù3ú V ü   õ(õ>ô  ÿ ÷ 6BKs- 'ûô Z ÷ -3 ÿ ùGúT÷Mù-ü  bôDù>ù-ÿ  ô   õ>ô A? ôú>ûôú>ûõ>ôWôù-üýþ÷ø j ÿ  9ù 6 I÷õTþÿ  `ú>÷ú>ûôúTûõ>ôWôô  ô 3 Lú>ÿ ÷  Põ>ô  ÿ ÷ 6 9ùjù-÷Nþô ? ôDþ u- ûôù-ô x üô  bôDù ÷;ö(øUôWù-û9ôWù Z ;õ>ô2÷ýú  ;ÿ  ôDþ]ý - ]õ>ô A? ôWøNô  Rúü9ù-ÿ  Pú>ûô noCp!py óóõ>÷ 66 Tû # 9þ]ú>ûô "3d_$[7bW]  ÷õ>ÿ ú>ûø  ùUþôDù 1 Iõ>ÿ ýrôDþ ô  ;õ  ÿ ôWõ  'û9ôWù-ô]õ>ôDù-ü  úTùUô Az ô  bú>ÿ  ô   vþôWøU÷  9ù3ú>õ A Lú>ôXú>ûô]ürù-ôX÷öú>ûô 9 õ  óûvý  ù-ôDþþ  ;ú  ù3ú>õ>ü  bú>üõ>ô#öF÷õ Z[7bDWC = UúTûô#ù-ô  I÷  9þUó  õ-ú÷;öú>ûô`ù3ú>ü9þ 7 'ôýõ>ÿ ô { Uõ>ôIó÷õ-ú÷  Uú>ûô#ýrôIû  ÿ ÷õ(÷öú>ûô ZnoCp!p ù3ú  Lú>ÿ ù3ú>ÿ  WùjöF÷õ ú>ûô q[7bDWy ;óó9õ>÷ 66 Tû <uc ù  ÷ú>ôDþÿ FU jõ>÷ór÷ù-ÿ ú>ÿ ÷ *: Xô  ;õ 1 ÿ ôIõ  9ú>û9ô jnoCp!p ÿ ù ! I÷øUóüú>ôWþdþÿ õ>ô  bú 1  ]öFõ>÷ø ú>ûôù-üý  õ  óûUÿ  P÷üõ V õ  ;ó9ûNý  ù-ôDþPþ  ;ú  ù3ú>õ>ü  bú>üõ>ô 6 i ôõ>ôWór÷õ-ú(õ>ôDù-ü  úTùGürù-ÿ  `ú  ÷øUôIú>õ>ÿ  Wù 6 'ûô >? 9õTù3ú øNôIú>õ>ÿ  `ÿ ù'ú>ûô C üøMýrôWõ÷ö O þþÿ ú>ÿ ÷  ú>õ>ÿ 36 ôDù'õTô A? ôDþ q û9ô  aõ>ô ? ÿ j|A 'û9ô2ù-ô  b÷ 6 9þaøUôbú>õ>ÿ  `ÿ ù'ú>ûô ø 3 ÿ øMüø } ô  ;ú>û4÷;ö nVoCp!p~ ù÷öjúTûô  ôWÿ  ûRýr÷õTù÷;ö |> ÿ ú>û9÷üú 2 b÷ü  ú>ÿ  wú>õ>ÿ 36 ô 2| i ô  ÿ   Gõ>ôIöôWõ ú>÷MúTûôDù-ôøUôbú>õ>ÿ  Wù  ù >y!3 9þ N]" Mõ>ôDù-órô  bú>ÿ  ô  6 E ';V1*7A q!6Z7-6 6Bj'1-6C13_ % 2-3<-(3 A(-  6  66 u!1' _ =24 ÿ  üõ>ô D7 IúTûô(øUô 323 9þ2ù3ú  9þ  õTþþô  Rÿ  Lú>ÿ ÷ 6 ÷;ö yV 9þ C¡" õ>ô T ÿ  ô < ûôõ>ôWù-ü  úTù ù-û÷ 3 4ú>û  Lú ¢6 ÷ 69% ôIõ  ô#úTô  rþùú>÷ #7E 'û9ÿ ùøNô  9ùú>û  Lúú>û9ô #A$ ù '& ô  ôIú>÷ 6q õ  ;óû ( ùóõ>ôDù-ô  Rú>ôDþaÿ  ù-ô  bú>ÿ ÷ "7  ÿ ù÷ N ôIõ A3 ôú>õ % ôWõTù-ôDþNöF÷õ _2 ÷þôDùórôWõô  ôWøUô  ú  ûôDù-ô`õ>ôDù-ü  úTù  ù  ô  < ùú>û÷ù-ôÿ  Xú>ûô`ó9õ>ô  ÿ ÷ü9ù ô A órôWõ>ÿ øUô  úTù#ù>û÷ % ú>û  ;úú>ûô  I÷ù-ú#÷;ö(órôWõ-ö÷õ>øUÿ  Pú>û9ô 2 ÷ 7 õ>ô A? ôWøNô  Rú#órôWõ#ô  ôWøUô  ú  ÿ ú>ûú>ûô j"d>$ î £A¤¥¦g§ ¨A©ªB« ¬¢A®A®A® á A£A©¯Z°q§ ± ²¥Z³ vä £A¯%´ ë ò 6ªgµ Râ ¶·A¸¸D¹º6»¼¹¸½'¾'»¿C½'À Á%»<ÂuºÃºà A®®A®ÄÅ3Å7Æ ÇÈÉ Ê ¾8½8ËÌA¾8½ÍZ¹ÎÏ ºÃ>ÐÑÒVÓAÔÕ Ö× ÐØ Ù A.3. Communications in Numerical Methods in Engineering 245 ÚOÛ<ÜÝEÞ(ßEÜTàáEâ(â<Üã7Ü*àãÛuäåEãäÛ<áuà>æçÛ(àèOáé7áãçê(ßEÜTàáEâ*Û<áæë/êáEìáuê<ãNÜéÚTçOÛuë/ãÞìCàîí ï ïð ñ ò òð ñ ó óAð ñ ô ôð ñ õ õð ñ ö ògï ÷ ó ô õ ñ ø ù ú-û ü8ýEþ6ü8ý1ü8û ÿ ü   ð  ò   óAð     ò  ó ï ï1ð ñ ò òð ñ ó óAð ñ ô ôð ñ õ õð ñ ö ògï ÷ ï ïð ñ ò òð ñ ó óAð ñ ô   8ý       ò  ó   "!$#&%('*)+'  ,-' ./,10 23541)7698:)36&)3  )<;=/,?>@# A&BCEDGF HIGJ*K L*MONPD$Q@QOJ*IDR/MOSTUDWVOX&SYZDGR7R3S7Q&L/DG[OF S]\$D$F ^OS_DT`L*MOS]a&^ONb[cS7J9IGd=J*S3V@aOS7N`S7aL9S7F S7N`S7aL/T K a@R3J*SDGTST7ef?LgD$F TIUT^OHHGST L/ThL*M@D(LgD`Y&S7NiD$a@YOK aOHUJ*S3V@a@S3N`S7aLgR3J*K L*S7J*K IGajR7D$ak[cS=^@TSYmlTK acR+S<L*MOS7J*S<n5K F F [cS=N`I&Y&S7J/D(L/S=DGYOYOK L*K IGacD$FoaOS7K HGM&[cIGJpT^O[Y&K \qK TK IGacr+e s MOK TpJ*ST^OF L5dIGJpL*M@Sbt1uwv=vxD$Q@QOJ*IDR/Mk^cTK aOH`L*M@S7TSwN`S3L*J*K R7TpK TgI$dyR3IGa@TK Y&S7J/DG[OF SwK aqL*S3J*ST L"TK a@R3S3L*MOS ^&L*K F K L z{I$d<L*MOK TULzqQcS]I$d<F IGa@HGST LiS7YOHGS|[OK TS7R+L*K IGaZT*R/MOS7N`S]D$a@Y{L*M@S3K J9S3}iR3K S7aLi^@TS]I$d<N`S7N`IGJ*z~McDGT QOJ*S7\qK IG^@TF zW[cS7S7a^OST L*K IGa@S7YeyO^OJL/MOS7JiYOS+L/DGK F TUn5K F Fp[cS|QOJ*STS3aqL*SY~K aN`IJ*SjYOS3QOL*MK aD_dIF F I(npK a@H T L*^@Y&ze  eg5"=h"A&f*""A  SkM@D\SiJ*S3\&K S7nS7YWTS7\GS7J/D$FyJ*S3V@aOS7N`S7aLbK Y&SDGTwK a~[cI$L/ML nhImD$acYmL*MOJ*S7SkY&K NUS7a@TK Ia@TbD$a@YWS+X&L*S7a@Y&SY IG^OJUQOJ*S7\&K I^@TbnhIJ*mIaTGS7F S+L*Ia&-[cDGTSYJ*S3V@aOS7N`S7aL`DGF HIGJ*K L*MONiT7e:OIGJbL*MOS]Gp A&B5CD$F HIGJ*K L*MON|DGa K N`QOJ*I$\GSYZD$acY\GS7J*za@D$L*^OJ/D$F"S7YOHGS3?Y@D(L/D~T L*J*^cR/L*^@J*S][@DTSYIaZHGJ/DGQOM@T9M@DTi[cS7S7aK aqL*J*I&Y&^cR3S7YZDGa@Y T L*^@Y&K SYoecAqS7\GS7J/D$FoQOJ*IGQcS7JL*K STI$dL*M@SwT^OQOQcIJL*K aOHiYOD(L/D`T L*J*^cR/L*^@J*SwD$J*S=Y&ST*R+J*K [cSYo@TMOI$n5K aOHbL*MOK TL*Ii[cS S+}9R3K S3aqL5K ajT L*IJ/D$HS"DGa@Y9L*K N`SwR3IT L(e s MOS<S3X&L*S3acTK IakIGdL*MOS<HJ/D$QOMk[@DTSYkYOD$L/DUT L*J*^@R+L*^OJ*S"dIJL*MOS= R3DTSbK L=K T<DGF T*IiQOJ/S7TS7aL/S7YcQ@J*IGQcIqTK aOHkDkR3IGJ*J*STQcIa@Y&K aOHidDR+S3?[@DGTSY]Y@D(L/DkT L*J*^@R+L*^OJ*Se  S`Y&K T*R3^@T*TgK a Y&S3L/D$K F$L*M@S1HJ/D$QOMcTcdIGJG{J*S3V@aOS7N`S7aLI$d&L*J*K DGaOH^OF D$L*K IGa@TDGa@Y"dIGJGUe7K a@D$F F z7aq^ON`S7J*K R3DGFGS3X&QS3J*K N`S7aL/T L*IUK F F ^cT L*J/D(L*S<L*M@S=K YOS7DTQOJ/S7TS7aL/S7YkK a9L*MOSwQ@D$QcS7J"D$J*S"QOJ*STS7aL*SYk^@TK aOHUL*MOSbGp AOBCD$F HGIJ*K L*M@Nje ÜåuèTêç o #é-áEâTÚáEìáEêãà >q ,g /,?+)3 .>|>&)7,g&//6|,?q&27?' +j 6|q)3?'"(]4h2 * 6q2ib8:)36&)+ )7,/o7)76('<h;<&*p*/¡¡¡¢$/£(¤$$ ¥p¦§m¨%q8y826('?*)3.*'h2#©:ª7£(«7¬¬ª$@)76qb$b¥p¦¥826$'?)7.*'12#h®h8¯-8®y¯-¡(«/¯-°°7°(±$# å c²+³3´µ¶ ·+¸3¹9º »½¼+¾+¾+¾5¿ ²+¸À"k¶ Á Â?´"à Fà 3²+À(ÄÅ %é ¹ÆqÇ È@É+ÊhÊË7ÌÍÎËÊ5ÏÐÍÑwÏÒ Ó(ÍÔoÌ7Õ3Ì7Õ ¼+¾/¾+¾ÖG×$×&Ø Ù?Ú7Û ÜÐÏÝÞ+ÐÏ?ß"Ëà-á Ì7Õpâ ã3ä1å+æ7ç èé â ê ë 246 Ap´endice A. Publicaciones ì íî(ïOîGð7ñ|ò ó:ôõoöO÷$ø î ùî$ú ó:ôõ@ûmó:üý]ó îïþ óyö$ó ôõ ù õôõü ú õ ð ÿ î  bí(î ÷ ï   Uí(î ó      ó:ý     !#" %$ ý !& '( )*'  + , !-* î ./02143145768:9;14<=90 >?@ABC/DE1 ÿ FF ÿ GIHKJIL ÿ M (ÿ N î O îoú CP ø î ùoî RQS?8TI60 U0 V ?/IUWCXC;V @AYZD)9/I9;[U0 V ?/\U)@UTI02U0 V ?/]U/I@-A?W 60 V ?/%A0 ;[U0292DV 9A1  P ^ _ù ' * ÷ ÿ FFN î $ îoú CP ø î ùoî ó !-* =_!-)`*' !&a" 7_ - !&)'! +b  * *î R31c?deQS?8TZ1<=90 >1SfCTTK1<=9g>1 BC/DE1 ÿ F)NhGZi  ÿ  L F $ M $ÿ jk î l îoú CP ø î ùî ÷ ð E !m`n Uî ÷ o :p î úî ó ' **"Ka*[ b ' b *"  Oq ý r$ q ý  +  s)!-* &_!&) î ./0[1431568&9;14<=90 >?@A7BC/DE1 ÿ FttGIuvKL Ohj)NM)OhOO î k î ø ** í(î þcî ÷  b 'w  î o î xI?T?W ?DV gUWXC;[UT>yxK>9?;zE1 Gíî o   Pe^ Wð )* ÷ ÿ FtN î h î {P gùî RXC;[UT>yxK>9?;zE1 ó  * o * P ÷ ÿ F)NO î N î p7[ @þcî/ñ *  ' + !-!& " 4*  7c*[ b ' b R"  +  P )* b " '* *î ZQS?8TK14Xc9?8-1 xK>9?;zmU/I@&fCT)TIW V gU0 V ?/A ÿ FFFGKJHKL hkMEFj î t îcþ K  5ïOî ý * b ' b *:'' + *" R +  s)c_ c !-)S!& E* î 4QS?8TI60 V /D ÿ FFkGE|I|ZL $ OkM $ kl î F îcï R } ó î ó *[P! +  ',I sE C"Z `s)"Z `~['' *c"K c +   'Z)  *C -*!&* ! +   + [  * î Zf75.[fcBR`7x(f75)2)E) ÿ FFFGO $$ MEOl)j î ÿ j îcï R } ó î ÷ ú CP ø î ùî(þ KE'4_!-)7"R* ! +   ' Z *,*%r a*w  î KfT)TK1S5768&1<=U0 >1 OjjjG HZuZL ÿ FkM)O ÿ t î ÿ+ÿ îoð b ò } í(î ï@î ÷ ú CP ø î ùoî ÷ ï  } ó î ÷ ï S ò :n Uî ó î ø  + &,*a7*[ b ' b *" *w :,*&_!-)    !&* î xZ.QKf7<CB`R7xyA9;V 9A Ojj ÿ GZKJ)J î ÿ O î ô  sn Uî úî ne* _!-)R,*aa c  }a, *' &"* ! +   '* *î .[f7<315768:9;1f/IUW 1 ÿ FtlG uKL hjlMEh ÿ $ î ú R + P  ) '  Ojjj í   o   Pr^ ~ð )* ÷ þ Z qî QS?8a8:6E/Z1568&9;1<=90 >1cBC/ED)/ED OjjjGI4L ÿ MEh C;[9[TU;[9@r6AV /ED-)E    Bibliograf´ıa [1] A.V. Aho, J. E. Hopcroft and J.D. Ullman,Data structures and algorithms, Addison Wesley, 1983. [2] D.N. Arnold, A. Mukherjee and L. Pouly,Locally adapted tetrahedral meshes using bisection, SIAM J. Sci. Comput. 22, 2 (2000), pp. 431–448. [3] I. Babuˇ ska and W.C. Rheinboldt,Error estimates for adaptive finite element computations, SIAM J. Numer. Anal., 15 (1978), pp. 736– 754. [4] I. Babuˇ ska and W.C. Rheinboldt,A posteriori error analysis of finite element solutions for one-dimensional problems, SIAM J. Numer. Anal., 18 (1981), pp. 565–589. [5] R.E. Bank and A.H. Sherman,The use of adaptive grid refinement for badly behaved elliptic partial differential equations, in Advances in Computer Methods for Partial Differential Equations III, R. Vichnevetsky and R.S. Stepleman, eds, IMACS, 1979, pp. 33–39. [6] R.E. Bank and A.H. Sherman and H. Weiser,Refinement algorithms and data structures for regular local mesh refinement, R. Stepleman et al., Scientific Computing, North-Holland, 1983, pp. 3–17. [7] E. B¨ ansch,Local mesh refinement in 2 and 3 dimensions, IMPACT Com. Sci. Eng., 3 (1991), pp. 181–191. [8] B.G. Baumgart,A polyedron representation for computer vision, Proc. AFIPS Natl. Comput. Conf., 44 (1975), pp. 589–596. [9] M. Berger,Geometry, Springer-Verlag, 1987. 247 248 BIBLIOGRAF´ IA [10] J. Bey,Tetrahedral grid refinement, Computing, 55 (1995), pp. 355–378. [11] J. Bonet and J. Peraire,An alternate digital tree algorithm for geometric searching and intersection problems, Int. J. Num. Meth. Eng., 31 (1991), pp. 1–17. [12] G.F. Carey,A mesh refinement scheme for finite element computations, Com. Math. App. Mech. Eng., 7, 1 (1976), pp. 93–105. [13] G.F. Carey, M. Sharma and K.C. Wang,A class of data structures for 2-d and 3-d adaptive mesh refinement, Int. J. Num. Meth. in Eng., 26 (1988), pp. 2607-2622. [14] P.G. Ciarlet,The finite element method for elliptic problems, in Studies in mathematics and its applications, J. Lion et al., Eds., NorthHolland, 1978. [15] P. Cignoni, C. Montani and R. Scopigno,A comparison of mesh simplification algorithms, Comp. & Graph., 22, 1 (1998), pp. 37-54. [16] T.J. Cross,A two-dimensional triangular mesh generator using the advanced front method, Dept. Civil Engineering, U. Swansea U.K., Technical Report CR/662/91. [17] L. De Floriani, P. Magillo and E. Puppo,Applications of computational geometry to geographic information systems, in Handbook of Computational Geometry, J.-R. Sack and J. Urrutia, Eds., Elsevier, 2000, pp. 333-388. [18] L. De Floriani, B. Falcidieno, G. Nagy and C. Pienovi,A hierarchical structure for surface approximation, Comp. & Graph., 8, 2 (1984), pp. 183–193. [19] H. Edelsbrunner,Algorithms in Combinatorial Geometry, SpringerVerlag, 1987. [20] J.M. Escobar,Generaci´on de mallas tridimensionales mediante triangulaci´on de Delaunay, Tesis Doctoral, Universidad de Las Palmas de Gran Canaria, 1995. BIBLIOGRAF´ IA 249 [21] L. Ferragut, NEPTUNO, Un c´odigo para el m´etodo de elementos finitos adaptativo, (en FORTRAN) Dept. Mat. Aplic. y Met. Inf., ETSI Minas, Madrid, 1987. [22] L. Ferragut, R. Montenegro and A. Plaza,Efficient refinement/derefinement algorithm of nested meshes to solve evolution problems, Comm. Num. Meth. Eng., 10 (1994), pp. 403–412. [23] L. Ferragut, R. Montenegro, G. Winter and A. Nu˜ nez,Accurate extraction of interconnect capacitances by adaptative mixed finite element method, Proceedings of the EUROMICRO 91, North–Holland (1991), pp. 61–68. [24] S. Fortune,Voronoi Diagrams and Delaunay Triangulations, in Computing in Euclidean Geometry, pp. 193–233, World Cientific, Singapore, 1992. [25] W.H. Frey, and D.A. Field,Mesh Relaxation: A new technique for improving triangulations, Int. J. Num. Meth. in Eng., 31 (1991), pp. 11211133. [26] J.P.S.R. Gago, D.W. Kelly, and O.C. Zienkiewicz,A posteriori error analysis and adaptive processes in the finite element method. Part II: Adaptive mesh refinement, Int. J. Num. Meth. Eng., 19 (1983) pp. 1621–1656. [27] H. George,Automatic mesh generation, J. Wiley, 1991. [28] M.T. Goodrich and K. Ramaiyer ,Geometric data structures , in Handbook of Computational Geometry, J.-R. Sack and J. Urrutia, Eds., Elsevier, 2000, pp. 463-489. [29] J.L. Gross, T.W. Tucker,Topological Graph Theory, John Wiley and Sons, 1987. [30] L. Guibas, J. Stolfi,Primitives for the manipulation of general subdivisions and the coomputation of Voronoi Diagrams, ACM Transaction on Graphics, 4(3) (1985) pp. 74–123. 250 BIBLIOGRAF´ IA [31] P.L. George, F. Hecht and E. Saltel,Automatic mesh generator with specified Boundary, Comp. Meth. in Mech. and Eng., 92 (1991) pp. 259–268. [32] R.L. Graham, D.E. Knuth, and O. Patashnik,Concrete Mathematics, John Wiley and Sons, 1989. [33] W. Hackbush and V. Trottengurg (eds),Multigrid Methods II, Lectures Notes in Mathematics, Springer–Verlag, 1982, Berlin. [34] W. Hackbush and V. Trottengurg (eds),Multigrid Methods II, Lectures Notes in Mathematics, Springer–Verlag, 1986, Berlin. [35] J. Hartmann and J. Wernecke,The VRML 2.0 Handbook, AddisonWesley, 1998. [36] D.M. Hawken, P. Townsend and F. Webster,The use of dynamic data structures in finite element applications, Int. J. Num. Meth. Eng., 33 (1992), pp. 1795–1811. [37] C.M. Hoffmann,How to construct the skeleton of CSG objects, in 4th IMA Conference on Mathematics of Surfaces, Bath, England, (1990). [38] C.M. Hoffmann,Geometric and solid modeling, Morgan Kaufmann Publishers, Inc., 1989. [39] H. Hoppe,Efficient implementation of progressive meshes, Comp. & Graph., 22, 1 (1998), pp. 27–36. [40] R. Hoppe,Bemerkung der Redaktion, Archiv Math. Physik (Grunert), 56, pp. 307-312, 1874. [41] H. Jin and R. I. Tanner,Generation of unstructured tetrahedral meshes by the advancing front technique, Int. J. Num. Meth. Eng., 36 (1993), pp. 1805–1823. [42] B. Johnston, J. Sullivan and A. Kwasnike,Automatic conversion of triangular finite element meshes to quadrilateral elements, Int. J. Numer. Meth. Eng., 31 (1991), pp. 67–84. [43] B. Kearfott,A proof of convergence and an error bound for the method of bisection in IRn, Math. Comp., 32 (1978), pp. 1147–1153. BIBLIOGRAF´ IA 251 [44] L. Kettner,Using generic programming for designing a data structure for poliedral surfaces, Comp. Geometry, 13 (1999), pp. 65–90. [45] I. Kossaczk´ y,A recursive approach to local mesh refinement in two and three dimensions, J. Comp. App. Math., 55 (1994), pp. 275–288. [46] P. Lindstrom, D. Koller, W. Ribarsky, L. Hodges, N. Faust and G. Turner,Real-time, continuous level of detail rendering of height fields, Procceedings of SIGGRAPH’96, (1996). [47] A. Liu and B. Joe,On the shape of tetrahedra from bisection, Math. Comp., 63 (1994), pp. 141–154. [48] A. Liu and B. Joe,Quality local refinement of tetrahedral meshes based on bisection, SIAM J. Sci. Comput., 16 (1995), pp. 1269–1291. [49] S.H. Lo,A new mesh generation scheme for arbitrary planar domains, Int. J. Num. Meth. Eng., 21 (1985), pp. 1403–1426. [50] S.H. Lo,Perspective projection of non–convex polyhedra, Int. J. Num. Meth. Eng., 26 (1988), pp. 1485–1506. [51] R. Lohner and P. Parikh,Three–dimensional grid generation by the advancing front method, Int. J. Num. Meth. Fluids, 8 (1988), pp. 1135– 1149. [52] The Mathworks Inc.,Partial Differential Equations (PDE) Toolbox User’s Guide , Matlab Online manuals. [53] J.M. Maubach,Local bisection refinement for n-simplicial grids generated by reflection, SIAM J. Sci. Stat. Comp., 16 (1995), pp. 210–227. [54] W.F. Mitchell,Optimal multilevel iterative methods for adaptive grids, SIAM J. Sci. Statist. Comput. 13 (1992), pp. 146–167. [55] M. Mantyla,An introduction to Solid Modeling, Computer Press, 1988. [56] J.R. Munkres,Elementary differential topology, Princeton Univ. Press, second edition, 1966. 252 BIBLIOGRAF´ IA [57] A. Mukherjee,An adaptive finite element code for elliptic boundary value problems in three dimensions with applications in numerical relativity, PhD. Thesis, Penn. State University, University Park, PA 16802, 1996. [58] D.E. Muller and F.P. Preparata,Finding the intersection of two complex polyedra, Theoret. Comput. Sci., 7 (1978), pp. 217-236. [59] E. Muthukrishnan, P.S. Shiakolas, R.V. Nambiar, and K.L. Lawrence,Simple algorithm for adaptive refinement of threedimensional finite element tetrahedral meshes, AIAA Journal, 33 (1995), pp. 928–932. [60] K. Nakahashi and D. Sharov,Direct surface triangulation using the advancing front method, AIAA–95–1686–CP, (1995). [61] J. Nievergelt and P. Widmayer,Spatial data structures: concepts and design choices, in Handbook of Computational Geometry, J.-R. Sack and J. Urrutia, Eds., Elsevier, 2000, pp. 725-764. [62] J. Pach (Ed.), New trends in discrete and computational geometry, Springer-Verlag, 1993. [63] M.A. Padr´ on,Un algoritmo de desrefinamiento en dimensi´on tres para mallas encajadas de tetraedros basado en el esqueleto, Tesis Doctoral, Universidad de Las Palmas de Gran Canaria, 1999. [64] R. Pajarola,Large scale terrain visualization using the restricted quadtree triangulation, in Proc. IEEE Visualization’98, 1998, pp. 19–26 & p. 515. [65] J. Peraire, M. Vahdati, K. Morgan and O.C. Zienkiewicz, Adaptive remeshing for compressible flow computations, J. Comp. Phys., 72 (1987), pp. 449–466. [66] A. Plaza,Algoritmos de desrefinamiento en mallados estructurados bidimensionales, Tesis Doctoral, Universidad de Las Palmas de Gran Canaria, 1993. [67] A. Plaza,Sobre el n´umero de tri´angulos generados por la partici´on 4T-LE, en preparaci´on.