Development of an integrated system for geometry modelling, mesh generation and data management for the analysis by the finite element method

Resumen

En este trabajo se describe el desarrollo e implementación de un sistema informático para el tratamiento de toda la información necesaria para un análisis por el Método de los Elementos Finitos o por otros métodos numéricos (diferencias finitas, volúmenes finitos, métodos de contorno, métodos de puntos, etc.).

Algunas de sus partes se refieren principalmente al diseño y organización de un sistema de estas características. En otras, se describen los nuevos algoritmos que ha sido necesario desarrollar para cumplir los objetivos propuestos.

Las diversas disciplinas que se describen a lo largo de la tesis se pueden diferenciar en:

• Organización del sistema. Consiste en definir el tratamiento de todos los datos de un análisis genérico de manera uniforme. También se dan criterios sobre el ordenamiento interno de los datos.
• Modelación geométrica. Se presenta una serie de algoritmos que se han desarrollado para tratar y modificar la geometría del modelo.
• Generación de malla. Se describen diferentes técnicas y algoritmos para poder realizar la generación de la malla.
• Adaptabilidad del sistema a diferentes análisis. Se describe como se ha diseñado y como se realiza la adaptación del sistema a un código de análisis cualquiera.

La implementación del conjunto de criterios y algoritmos descritos a lo largo de esta tesis, ha permitido la creación de un sistema que da soporte al proceso de análisis mediante métodos numéricos de modelos, tanto a nivel académico como industrial.

1 Introducción

1.1 MOTIVACIÓN

Con el desarrollo de los ordenadores en los últimos años, se ha llegado a un momento en el que es posible realizar análisis de modelos extremadamente complejos. Además, los programas de cálculo y simulación por métodos numéricos han ido sofisticándose con el paso del tiempo para modelizar de manera más aproximada los fenómenos reales. Ello conlleva que la complejidad de los datos que requieren haya ido aumentando. Otro efecto de esta evolución radica en que las mallas necesarias deben cumplir unas especificaciones muy determinadas y ser capaces de modelar situaciones muy especiales de la geometría¹ en el contexto del problema bajo estudio.

El resultado de esta necesidad ha sido que, actualmente, el proceso que generalmente consume más tiempo-hombre en el conjunto de un análisis es la preparación de datos y la generación de la malla.

El segundo problema reside en que en muchos casos ni tan solo es posible llegar a obtener una malla y unos datos que cumplan los requerimientos de un análisis determinado.

Otra necesidad que aparece es la de conectar de manera adecuada el sistema de preproceso con el de análisis, con el objetivo de minimizar el tiempo y la dificultad de la transferencia de datos.

Finalmente, ha llegado el momento de que los métodos numéricos traspasen la barrera de la investigación y se introduzcan de pleno en la producción industrial. Un sistema sólo será usable industrialmente cuando se supere toda la problemática de la interacción con el usuario así como cuando una persona, con conocimientos técnicos adecuados, sea capaz de utilizar un programa sin excesivas dificultades prácticas.

(¹) No es suficiente con una malla que describa al sólido. Hay que pensar en otros factores como contactos, juntas, simetrías espaciales, multi-volúmenes y otros que requieren una malla especializada para realizar con éxito el análisis.

1.2 ANTECEDENTES

Los primeros estudios de modelación geométrica se remontan a los años 60-70 y se originaron en las grandes compañías fabricantes de automóviles. En un inicio, estas compañías fueron las creadoras y las usuarias, al mismo tiempo, de estos primeros programas de diseño asistido por ordenador (CAD). Más tarde, los programas de este tipo se han ido extendiendo a todos los ámbitos de la ingeniería y del procesamiento industrial hasta el extremo de que, hoy en día, prácticamente cualquier diseño o mecanización en la ingeniería se realiza mediante este tipo de programas. Puede considerarse que se ha llegado a un estado del arte, tanto en investigación como en programas comerciales, que cumple con eficacia su cometido.

A este respecto, deben tenerse en cuenta los trabajos de Bézier, Coon, Gerald Farin [1], Foley, van Dam,Feiner y Hughes [2] y otros que han permitido que, en unos pocos años, todo el conocimiento matemático sobre geometría computacional se haya extendido y sea ampliamente usado en todos los programas de modelación geométrica tridimensional.

La generación de malla ha sido durante los últimos 15 años y sigue siendo hoy en día objeto de investigación por parte de multitud de grupos tanto en la universidad como en el ámbito de la explotación comercial. Actualmente existen ya metodologías que funcionan razonablemente bien para algunos de los casos en que se subdivide esta disciplina. Otros apartados, como generación no estructurada de hexaedros, capas límite para análisis de fluidos viscosos gobernados por la ecuación de Navier-Stokes y otros pueden considerarse como no resueltos o, al menos, como muy mejorables. Además, el proceso de relación entre geometría y generación de malla provoca una problemática que requiere procesos de adaptación del modelo geométrico, corrección de imperfecciones asociadas y preparación para la generación, que puede considerarse como un problema no resuelto de manera totalmente satisfactoria.

Los primeros trabajos sobre generación de malla no estructurada bi y tridimensional se deben a R. Lohner [3], J. Peraire [4] y J. Peiró [5]. Más tarde, el número de investigadores y grupos que se han dedicado a las mejoras en las metodologías de generación ha ido aumentando progresivamente. Entre ellos pueden citarse a P.L. George y H. Borouchaki [6], Buscaglia, Gustavo C. y Enzo A. Dari [7], Peter Möller y Peter Hansbo [8], R. Kreiner [9], Jas Frykestig [10], A. Rassineux, J-M. Savignat, O. Stab y P. Villon [11] y muchos otros.

Los programas comerciales que permiten generación de malla comienzan a ser capaces de tratar modelos geométricos complejos con poca intervención del usuario. Sin embargo, no pueden considerarse aún del nivel de calidad que atribuíamos a los programas de CAD del apartado anterior. Adolecen de una cantidad demasiado alta de problemas y dificultad de uso. Por otro lado, son utilizables para problemas muy determinados y especializados y por tanto, no aptos para multitud de análisis. Los problemas referidos anteriormente se centran en la dificultad de obtener una malla que cumpla con unos requisitos prefijados por el usuario. Ejemplos de programas comerciales de generación de malla son: PATRAN [12], HyperMesh [13], Femgen/Femview [14], FEMAP [15], ICEM CFD [16] y otros.

En relación a la adaptabilidad, entendiéndose como la capacidad de relación entre un sistema de preproceso y un programa de cálculo, la mayoría de los programas de preproceso y postproceso existentes en el mercado no permiten su adaptación fácil a códigos determinados de análisis. Si se acepta que deben existir gran cantidad de códigos de cálculo, tanto para poder realizar investigación con ellos como para resolver problemas especializados, y se reconoce que no es posible que cada código de análisis contenga su propio sistema de pre-postproceso, se puede concluir que es necesaria la creación de un sistema de pre-postproceso que tenga gran capacidad de adaptación para dar soporte a los códigos antes mencionados.

1.3 OBJETIVOS

A partir de las motivaciones y criterios ya mencionados, el objetivo de la tesis es la creación de un sistema de preproceso con capacidad de tratamiento de modelos geométricos complejos, que permita la generación de mallas que se adapten a todas las tipologías de análisis y con posibilidad de tratamiento de todos los datos asociados a dichos análisis.

El sistema debe ser capaz de crear un entorno de trabajo que permita a un usuario preparar todos los datos que requiera un cálculo y minimizar el coste de tiempo que requiere esta preparación. Asimismo, la curva de aprendizaje del sistema debe resultar en un tiempo mínimo de iniciación y en una alta eficiencia de proceso al final de este período inicial.

Finalmente, se debe proporcionar una alta adaptabilidad entre el preproceso y el programa de simulación numérica. Se considerará óptimo si este proceso de adaptación puede realizarse exteriormente a ambos programas y con un coste de tiempo que se mantenga dentro de límites razonables.

1.4 PRESENTACIÓN DE LOS CONTENIDOS DE LA TESIS

Esta tesis presenta un conjunto de materias cuyo hilo común reside en que han sido necesarias para la implementación de un sistema de preproceso para el análisis por métodos numéricos (elementos finitos, diferencias finitas, volúmenes finitos, etc.), cuyo nombre es GiD. Cada capítulo cubre una parte de este conjunto común, bastante diferenciada de las otras y que pretende ilustrar las técnicas y algoritmos que han sido necesarias para llegar a los objetivos propuestos.

En general, cada capítulo incluye una serie de definiciones y consideraciones generales. Seguidamente describe el conjunto de algoritmos implementados, intentando hacer hincapié en los que han aportado conceptos e ideas novedosas al tema propuesto. En muchos casos es necesario referise a conceptos e ideas que forman parte del estado del arte para permitir una comprensión global. Al final de cada capítulo se incluye un apartado donde se especifica que ideas, algoritmos y metodologías se consideran como aportación de esta tesis al tema descrito.

A continuación se incluye una breve descripción de cada capítulo y de la temática que aborda.

El Capítulo 2 de la tesis describe la organización general del sistema que se ha desarrollado como preproceso para el análisis mediante métodos numéricos. En él se introducen algunos conceptos importantes relativos a la programación y estructuras de datos que se han implementado. También se proponen criterios para la definición de la interfaz de usuario y se especifica el tratamiento del sistema de unidades.

El Capítulo 3 aborda el tratamiento de la modelización geométrica. Para su correcta comprensión es interesante consultar el apéndice dedicado a la descripción matemática de las NURBS. Dicho capítulo da, en una primera parte, un conjunto de definiciones y criterios sobre las diversas formas de tratar la geometría. Seguidamente introduce el tema de la importación de modelos geométricos y describe un conjunto de algoritmos que se han implementado para la corrección de la geometría y su adaptación a la generación de malla. Finalmente, describe un conjunto de algoritmos relativos al tratamiento geométrico, cuya característica común consiste en que han requerido de un cierto estudio e investigación en contraposición a otros algoritmos que ya pertenecen al estado del arte. Algunos de estos algoritmos se describen en el mencionado apéndice sobre las NURBS.

El Capítulo 4 trata el tema de la generación de malla. En su inicio, se dan algunas consideraciones de tipo general sobre en que consiste dicha generación y en las particularidades y casuísticas que contiene. Seguidamente se procede a la descripción de las técnicas usadas para la generación. Durante todo el capítulo, se intenta hacer hincapié en los algoritmos que han sido creados o mejorados en la presente tesis y que, en algunos casos se mezclan con algoritmos ya existentes y extraídos de la bibliografía. Al final del capítulo se incluye un apartado en el que se precisan las aportaciones de esta tesis a las técnicas y metodologías de generación de malla.

El Capítulo 5 se refiere a la adaptación del sistema para un tipo de análisis. En él se describe qué se considera por adaptación y como se ha tratado dentro del sistema implementado. Asimismo, se crea el paralelo entre la abstracción de datos necesaria para un sistema genérico de preproceso y su implementación práctica.

El Capítulo 6 describe un conjunto de ejemplos de aplicación. Su objetivo consiste en ir introduciendo una serie de modelos, que crearon en su momento la necesidad de implementar nuevos algoritmos que pudieran tratarlos correctamente. Al mismo tiempo, dan una idea de que tipologías y en que nivel de dificultad se pueden resolver con los algoritmos desarrollados e implementados en el sistema GiD, que es el resultado más notable de esta tesis.

El apéndice A versa sobre la descripción matemática de las NURBS. Este apéndice contiene la descripción de las principales entidades que se consideran hoy en día imprescindibles para un sistema de modelado tridimensional. Al mismo tiempo, se incluyen algunos de los algoritmos más ampliamente conocidos y necesarios para el tratamiento de dichas entidades geométricas. La inclusión de este apéndice proviene del hecho de que muchas de las definiciones y algoritmos contenidos en él son de necesario conocimiento para comprender correctamente el Capítulo 3 y algunas partes del 4.

2 Organización del sistema GiD

2.1 INTRODUCCIÓN

En este capítulo se van a describir las principales pautas de diseño y de organización empleadas en el sistema de preproceso desarrollado en la tesis (denominado en adelante GiD). Es necesario tener en cuenta que dada la complejidad y las posibilidades que este ofrece, un diseño correcto es básico para que pueda ser usado de manera eficiente. A esto hay que añadir el hecho de que debido a que es un tipo de sistema orientado hacia el usuario final, la coherencia de uso cobra especial importancia.

Se va a describir la organización a nivel del código interno en cuanto al empleo de herramientas y utilidades externas así como librerías gráficas o de otros tipos.

Otra de las partes que se tratarán en este capítulo será la interrelación con el usuario¹ y las técnicas para definir un uso coherente y una metodología sencilla de aprendizaje.

Se va a dividir el conjunto de datos necesarios para el análisis en una serie de subconjuntos lógicos y funcionales que se han implementado en el sistema GiD.

Se especificarán con detalles las estructuras de los datos perceptibles por el usuario y se describirán los criterios lógicos para la adopción de tales estructuras.

Además, se describirá la metodología para tratar dentro del sistema el tema de las unidades y su adaptación a códigos de análisis que pueden o no incluir el tratamiento de éstas.

(¹) Usualmente llamada interfaz de usuario o user interface.

2.2 LENGUAJE DE PROGRAMACIÓN Y LIBRERÍAS

El lenguaje principal de programación ha sido C++. Por ser un lenguaje orientado a objetos, facilita el desarrollo de códigos que no tengan un flujo lineal de ejecución. En comparación a un típico programa de cálculo por ordenador, donde hay una secuencia ordenada de operaciones a realizar, en este caso el flujo, es fuertemente no continuo. Esto se debe al hecho de que está guiado por la interacción con el usuario. La consecuencia de ello es que interesa más ordenar el código por unas entidades abstractas llamadas objetos y que representan conceptos reales en lugar de definir el habitual conjunto de funciones.

Para facilitar el trabajo de programación y evitar la repetición de trabajo ya existente, se han utilizado diversas librerías. Se entiende por librerías, un conjunto de funciones escritas en un lenguaje compatible con el del programa y que están especializadas en realizar un trabajo concreto y muy bien definido. Las librerías usadas han sido:

• OpenGL : Librería para dibujo en 3D. Dadas unas primitivas de dibujo como líneas, superficies o triángulos, y dados un tipo de visualización y unos parámetros de posicionamiento, dibuja las primitivas en pantalla.

• Tcl-Tk : Librería para el uso de scripts¹ y para la construcción de la interfaz de usuario. Sirve para crear todos los menús, botones y ventanas que interaccionan con el usuario. La posibilidad de ejecución de scripts permite que una tercera persona pueda añadir funcionalidad al programa para realizar tareas concretas o para crear nuevas ventanas.

• Gettext : Librería de soporte para internacionalización. Facilita la traducción del programa a otros idiomas.

• Otras: Se han usado diversas librerías y código existente para operaciones tales como: escritura de ficheros de impresión, creación de ficheros de vídeo, etc.

Por ser Tcl-Tk una librería que opera con scripts, se puede desarrollar todo el código que gobierna la creación de la interfaz de usuario mediante este nuevo lenguaje que no necesita de compilación y por tanto implica un desarrollo más rápido y sencillo. Por la misma razón de que la compilación no es necesaria, una tercera persona puede añadir nuevas opciones y ventanas al programa para especializarlo en una tarea determinada. Esta unión entre el código en C++ y el código en Tcl requiere de una serie de funciones de interrelación que permitan el trasvase de información de una parte a la otra.

Una de las mayores ventajas de estas librerías, es que están disponibles para diversos sistemas operativos. Ello implica que el programa puede funcionar correctamente en toda una variedad de plataformas tales como con X-windows, MS-Windows , etc.

(¹) Un script es una secuencia de comandos que realizan acciones y que no necesita ser compilada.

2.3 CRITERIOS PARA LA CREACIÓN DE LA INTERFAZ DE USUARIO

Figura 1: Interfaz de usuario centrada en una ventana única.

Figura 2: Interfaz de usuario con ventanas múltiples.

En un programa dirigido a un uso directo con el usuario, tan importante como las opciones que este tiene y los algoritmos en él implementados, es la interrelación con la persona que debe usarlo. Si consideramos como objetivo de uso la creación de todos los datos de un modelo para su posterior análisis en el mínimo tiempo y con la máxima comodidad posible¹, la comunicación persona-máquina cobra especial importancia.

A lo largo de los años se han ido definiendo una serie de estándares que marcan unas reglas sobre como deben ser los menús, ventanas, cuadros de diálogo y todos los elementos que interaccionan con el usuario. Estos estándares están fuertemente relacionados con el sistema operativo que los utiliza y no son intercambiables con otros sistemas operativos. Ejemplos de ello son:

• Sistema operativo Macintosh: El precursor en ordenadores personales. La apariencia de todos los programas que funcionan en él es muy similar y bien definida.

• Sistema operativo MS-Windows: A partir de los conceptos definidos por el anterior ha creado sus librerías MFC² y unas reglas también muy definidas de como deben ser todos estos elementos.

• X-Windows y Motif: Proporcionan también un estándar de creación de widgets³. En este caso las reglas de creación son un poco más libres.

Partiendo de la base de que se quiere crear un sistema no dependiente del sistema operativo y teniendo en cuenta todos estos factores y estándares, se han establecido una serie de reglas en las que se ha basado la posterior creación de la interfaz de usuario:

1. Debe ser posible trabajar con el programa tanto en un entorno de ventana única como en uno de múltiples ventanas.

Comentario: Según la capacidad del ordenador, tamaño de la pantalla, tipo de modelo de análisis o preferencias del usuario, puede ser más interesante tener una única ventana de trabajo con unos pocos menús y botones desde donde se pueda acceder a toda la funcionalidad del programa. En otros casos, el acceso a los comandos puede ser más intuitivo o cómodo mediante el uso de muchas ventanas. La preferencia por un esquema u otro debe dejarse al usuario.

2. Las ventanas en general no deben tener tamaño fijo y deben ser, hasta cierto punto, adaptables a todos los tamaños.

Comentario: Fijar el tamaño de las ventanas secundarias de entrada de datos asegura que estas tengan siempre una forma correcta. El grave inconveniente es que restringe la posibilidad de que puedan adaptarse a la disposición preferida por el usuario. Se acepta por tanto que pueda haber una cierta incorrección en la presentación a cambio de la mayor libertad de posicionamiento.

3. El modo de uso debe ser bastante similar a los estándares antes consignados.

Comentario: Puede sacarse un denominador común de las diversas formas de representación. Aprovechar estas formas de trabajo ya conocidas por el usuario reduce el tiempo de aprendizaje.

4. Toda la entrada de datos por parte del usuario, aunque pueda ser introducida de muchos modos (por ventanas, botones, menús, etc.), debe estar centralizada internamente en un solo lugar que pasa a ser el centro de flujo del sistema.

Comentario: De esta manera los diferentes estilos de entrada de datos producen los mismos resultados y pueden ser tratados de manera unificada. A partir de ello se facilitan una serie de operaciones tales como: creación de ficheros batch , utilidad de deshacer y otras. Ver [17] para más detalles sobre el tema.

5. Ser agradables de uso (estética y funcionalmente).

Comentario: Un usuario se pasa muchas horas utilizando un mismo programa. Es básico que le resulte adecuado su funcionamiento y disposición de sus partes.

6. Minimizar la curva de aprendizaje.

Comentario: Lo que resulta intuitivo para el usuario acostumbra a ser producto de un gran esfuerzo en la programación y diseño de un sistema.

7. Dar la máxima potencia y velocidad de uso a los usuarios expertos.

Comentario: En programas de uso no ocasional, como puede ser un sistema para análisis numérico, cobra especial peso la maximización de la eficiencia para el tratamiento de modelos complejos.

A partir de los criterios mencionados y de las posibilidades que ofrece la librería Tcl-Tk en cuanto a la creación de widgets, se han definido todos los componentes gráficos del sistema. Se puede comprobar que el mayor condicionante ha sido dar la libertad al usuario para que elija el sistema de trabajo que mejor se adapte a sus necesidades.

(¹) La calidad de estos datos, como podría ser la calidad de la malla generada, se discute en otros capítulos.

(²) Microsoft Foundation Class Library.

(³) Se llama widget a cualquier elemento gráfico de interacción con el usuario como menús, ventanas, botones, cuadros de diálogos, etc.

2.4 DESCRIPCIÓN BÁSICA DEL SISTEMA Y DE SU CONFIGURACIÓN

Dada la utilidad requerida al sistema: generar todos los datos para un análisis por ordenador, el programa debe tener conocimiento de los datos que requiere ese análisis en concreto. Pero al mismo tiempo el programa debe ser capaz de generar datos para cualquier análisis. La manera de resolver esta paradoja es que, en principio, el programa no tiene información sobre ningún dato particular del análisis. Sólo tiene las utilidades comunes a todo tipo de cálculos como son la modelación geométrica, generación de malla, y capacidad genérica para tratamiento de datos. En este momento, no hay nada que se pueda asociar a una fuerza aplicada o a una velocidad inicial. Pero existe la posibilidad de leer una determinada información de configuración diferente para cada código de análisis que le indica toda la información que necesita saber, y que en su momento el usuario deberá tratar.

El proceso se basa en que para cada código de análisis que se desee tratar, deben existir un conjunto de ficheros de configuración. Estos indican al programa qué datos debe introducir el usuario y cómo debe escribirse el fichero de datos que posteriormente leerá el código de análisis. En apartados posteriores se explicará con más detalle cuales son estos datos.

Una vez se tiene la configuración para diversos programas, se puede escoger en cualquier momento con qué análisis se va a trabajar de entre todos los instalados. Esto provocará que algunas de las ventanas y menús cambien para permitir al usuario introducir los datos que se necesiten para el cálculo.

2.5 AGRUPACIÓN Y ESTRUCTURA DE LOS DATOS

Observando los datos necesarios para un análisis cualquiera por métodos numéricos, se puede subdividir el total de la información en una serie de conjuntos funcionales que se describen a continuación:

• Malla que defina la geometría del modelo.

• Condiciones de contorno y, en general, toda la información asociada a los nodos o elementos de la malla.

• Materiales.

• Datos generales del problema.

• Datos generales de intervalo.

Estos subgrupos engloban toda la información necesaria y se van a describir en apartados posteriores.

2.5.1 Malla

En todo problema de elementos finitos, volúmenes finitos, diferencias finitas e incluso de estructuras de barras se necesita una malla que defina la forma geométrica del modelo. Esta malla puede estar compuesta de un solo tipo de elemento o de una mezcla de elementos de diversas tipologías. También puede contener elementos especiales como son los de junta o los de contacto. Los elementos de barra pueden considerarse como elementos 1D de dos nodos.

2.5.2 Condiciones

Se entiende por condiciones todos los datos que van asociados a entidades geométricas o de la malla. Entre ellos figuran las condiciones de contorno y las condiciones iniciales del problema. Pero en este subconjunto se puede incluir cualquier otro dato como podrían ser secciones o dimensiones de los elementos, masas ficticias de los nodos, etc.

2.5.3 Materiales

En su sentido más general, un material podría estar incluido en el conjunto condiciones porque cumple la definición de ser un grupo de datos asociado a los elementos. Pero para este caso concreto interesa diferenciarlo por la evidente significación e importancia que tienen los materiales en la mayoría de análisis y simulaciones. A nivel de tratamiento de datos, un material también es un conjunto de datos aplicado a los elementos, pero con la particularidad de que hay una base de datos de materiales con nombre, que pueden ser modificados, y la misma base de datos puede ser ampliada por el usuario.

2.5.4 Datos generales del problema

Se incluyen en este subgrupo todos los datos que son únicos para todo el problema. Ejemplos de ello podrían ser el tipo de análisis, solver¹ a usar, el tipo de integración para los elementos, etc.

(¹) Se llama solver al método de resolución de sistemas de ecuaciones.

2.5.5 Datos generales de intervalo

Son un conjunto de datos que pueden ser diferentes para cada uno de los intervalos. A continuación se incluye la definición de intervalo.

Definición de intervalo: Se entiende por intervalo un subconjunto de información que permite encapsular un conjunto de datos en grupos. Puede haber un número arbitrario de ellos en un modelo. Al iniciar un nuevo modelo habrá siempre un intervalo definido y el usuario puede crear más. Dentro de cada uno de estos intervalos, cambiará la información de condiciones y los datos generales de intervalo. El resto de la información es independiente de él.

Un caso típico de uso de los intervalos es cuando, dado un análisis transitorio, puede ser interesante aplicar al principio un tipo de condiciones de contorno (como podrían ser fuerzas aplicadas) y más adelante cambiarlas por otras. En este caso las condiciones serían diferentes para cada uno de los intervalos. También podría ser interesante para un problema del mismo tipo, cambiar la velocidad de avance en diferentes momentos del análisis. Esto se podría hacer creando los parámetros que controlan esta velocidad de avance como datos generales de intervalo.

2.6 SISTEMAS DE UNIDADES

Tradicionalmente, la mayoría de los programas de análisis por ordenador no han tenido en cuenta las unidades de las diferentes magnitudes y propiedades que se introducen para el cálculo. El único requerimiento ha sido, y sigue siendo en la mayoría de los casos, que todos estos datos estén introducidos en unas unidades compatibles entre ellas, sin necesidad de especificar cuales son en cada caso. Esto traslada la problemática y el control de errores sobre estas unidades hacia el usuario. Aunque es evidente que el usuario debe controlar perfectamente la coherencia de todos los datos, se obvia la necesidad de poner unos controles que faciliten la detección de errores.

Por otra parte, se puede comprobar fácilmente que el hecho de que no se realice control de unidades en el programa de análisis, es fuente de abundantes errores para el usuario. Incluso a nivel de visualización de resultados, la identificación por parte del programa de cada valor con su unidad correspondiente, facilita la comprensión de esos resultados y su interpretación posterior.

Se deduce de todo ello que es necesario que el pre y postprocesador de un programa de cálculo se ocupe de toda esta problemática y haga un control adecuado de errores y de la coherencia sobre los datos que introduce el usuario.

2.6.1 Problemática

La carencia de unidades en el programa de cálculo provoca el primer problema, ya que éste no tiene manera de indicar en el fichero de postproceso las unidades en que está expresado cada resultado. Deben buscarse medios adicionales de comunicación entre pre y postproceso para solventar esta dificultad. Hay que tener en cuenta la casuística de si el programa se ejecuta desde dentro del entorno (máximo control), o se escribe el fichero de datos, y se ejecuta el cálculo desde el exterior (mínimo control).

Otro problema es que para obtener la máxima funcionalidad en la elección de unidades distintas, es necesario que el preprocesador mantenga una tabla interna de unidades para diferentes magnitudes. Como la pretensión del programa es que se pueda adaptar a cualquier tipo de análisis, resulta evidente que puede haber algunas unidades o incluso magnitudes, que no estén contempladas en su interior. Por tanto, debe permitirse al tipo de problema, definir sus propias unidades y magnitudes que se sumarían a esta tabla general.

Finalmente, debe evitarse la confusión al usuario entre las unidades de trabajo, que pueden ser no compatibles entre ellas (siempre que esté claramente especificado), con las que se utilizan en el análisis, que sí deben serlo.

2.6.2 Esquema escogido

En el esquema considerado como más óptimo ¹, el preproceso tiene como información interna un conjunto de magnitudes, y diversas unidades para cada una de ellas. Para unidades compuestas, o sea que se forman como una combinación de las básicas, debe especificarse en primer lugar la combinación a partir de las básicas. En el caso de las unidades no compuestas, la unidad de referencia y la básica coinciden. La información se puede esquematizar según la Tabla 1.

Tabla. 1 Introducción de unidades y magnitudes físicas en el programa

Magnitudes		Reference unit	Basic unit	other units
MAGNITUDE:	length	m	m	100cm	…
MAGNITUDE:	mass	kg	kg	$10^{3}$gr	…
MAGNITUDE:	strength	${\frac {\mathrm {kg} \cdot {}\mathrm {m} }{\mathrm {s} ^{2}}}$	1N	$\approx {10}^{-1}{\textrm {kp}}$	…
MAGNITUDE:	pressure	${\frac {\mathrm {kg} \cdot {}\mathrm {m} }{\mathrm {m} ^{2}\mathrm {s} ^{2}}}$	1Pa	$\approx {10}^{-5}{\frac {\textrm {kp}}{{\textrm {cm}}^{2}}}$	…
MAGNITUDE:	temperature	d_Celsius	d_Celsius	${\frac {5}{9}}(\mathrm {d_{F}arenheit} -32)$	…
…

El programa de por si, ya tendrá un conjunto elevado de magnitudes y unidades, que podrá incrementarse para cada tipo de problema mediante el formato anteriormente descrito.

Por otra parte, cuando en el tipo de problema se definan propiedades de materiales que requieran unidades, se incluirá en la definición un indicador de que la propiedad es de este tipo y el valor del campo será el valor numérico más el nombre de la unidad, expresada según una de las unidades ya definidas. Por tanto, un ejemplo de una propiedad expresada según este criterio sería:

QUESTION: Young_modulus#UNITS#
VALUE: 2.1e11Pa

El indicador #UNITS# sirve para identificar a esta propiedad como del tipo “propiedad con unidades”, y así poder distinguir el valor 2.1e11Pa como formado por un valor numérico 2.1e11 y una unidad Pa.

(¹) Más óptimo tanto para el usuario, como para el que define el tipo de problema.

2.6.3 Control de unidades por el usuario

Únicamente con los datos del apartado anterior, ya sería posible trabajar con unidades dentro del programa. Con este esquema, las unidades de trabajo para cada propiedad y para la magnitud de longitud del modelo serían fijas. Si se supone que todas las propiedades de los materiales del tipo de problema vienen dadas en unidades compatibles y que la unidad usada en el modelo geométrico es también compatible con las anteriores, entonces el problema está bien definido.

El inconveniente de esta solución es su excesiva rigidez, ya que la unidad de longitud del modelo geométrico queda definida a priori sin que pueda ser cambiada y que las unidades de las propiedades también tienen que expresarse en unas unidades fijas. El primer inconveniente es el más grave, porque un modelo puede venir ya dado en unas unidades distintas y puede ser una dificultad adicional el escalarlo. Por otra parte, en general se trabaja con más comodidad en unas unidades tales que los números resultantes sean de un orden cercano a la unidad¹. El usuario está también acostumbrado a unas unidades determinadas, o tiene datos en esas unidades, y representa una dificultad adicional el cambio a las unidades exigidas.

Por todo ello, se hace necesario dar facilidades adicionales de control al usuario. La técnica propuesta consiste en un conjunto de preferencias donde el usuario indique las unidades en que desea trabajar. Este conjunto podría subdividirse en tres grupos diferenciados² que serían los siguientes:

1. En el primero de ellos se declararía la unidad de longitud en que está expresado el modelo geométrico. Esta unidad es básica porque es la especificación de como es el modelo, que en ningún caso va a escalarse.
2. En el segundo grupo se especifica que unidades se quieren utilizar para las propiedades de los materiales. Una de estas propiedades será también la longitud y su unidad no tiene porque coincidir con la del modelo.
3. El tercer grupo especifica en qué unidades van a enviarse los datos desde GiD hacia el módulo de cálculo. Una de sus unidades será también la de longitud y deberá estar fijada e igual a la unidad definida en el primer grupo.

La diferencia entre el primer y segundo grupo es grande porque en la especificación de la longitud del primer grupo, el usuario declara que la geometría del modelo está definida en esa unidad. Esta es una información que si se cambia a otra unidad que no sea la real, el resultado final será incorrecto porque el programa no puede comprobar este dato.

En cambio, las unidades del segundo grupo se utilizan solamente a efectos de edición y visualización. Como GiD conoce en que unidades están expresadas las propiedades originales, podrá escalarlas adecuadamente para que en GiD se visualicen en el sistema elegido. Incluso si el usuario cambia una propiedad en un momento dado y la introduce en las unidades activas en ese momento, si más tarde se cambian las unidades de visualización, la unidad se escalará correctamente de forma automática. Para evitar confusiones, es básico que toda propiedad tenga indicado visualmente en qué unidades está expresada y que también se indique muy claramente la unidad de la geometría del modelo.

La razón de introducir está complejidad adicional de que la unidad de longitud del modelo pueda ser diferente de la de visualización de propiedades radica en el hecho siguiente: el modelo normalmente está en unas unidades de trabajo no modificables, bien porque procede del exterior (importación de CAD o de otro modelo), o bien porque se trabaja de manera más apropiada con unas ciertas unidades. Por otra parte, las propiedades de los materiales se conocen en un sistema de unidades determinado y cambiarlas a otras no habituales puede ser una fuente de errores. Dos ejemplos típicos de ello serían:

• En algunas ramas de la ingeniería de procesos industriales, los modelos se trabajan típicamente en milímetros y el sistema de unidades para los materiales acostumbra a ser el Sistema Internacional, que tiene como unidad básica de longitud el metro. Modelar la geometría en metros sería muy poco práctico y fuera de costumbre. Por el otro lado, introducir los materiales con unas unidades que incluyeran al milímetro sería, como mínimo, poco ortodoxo además de proclive a errores.
• En ingeniería civil está muy extendido expresar las presiones en ${\textstyle {\frac {\mathrm {kp} }{\mathrm {cm^{2}} }}}$ y, según el tamaño de la estructura, puede ser de interés definir sus dimensiones en metros. En este caso también tendría interés el uso de ambas unidades.

El tercer grupo de unidades son las que se van a usar para enviar los datos al programa de cálculo. La razón de que la unidad de longitud esté fijada al mismo valor que la dada para definir el modelo, se basa en que es preferible evitar el escalado de los nodos de la malla. Este escalado, además de ser costoso en tiempo de computación, es propenso a provocar errores y confusiones. Es mucho más sencillo escalar las propiedades que tengan unidades. En general no importará en que unidades se envíen las propiedades al programa de cálculo, ya que en el postproceso habrá facilidades similares para volver a hacer el cambio a las unidades deseadas (ver Apartado 2.6.4). Se permite que el usuario tenga control sobre este tercer grupo de unidades con la intención de que, en casos especiales, pueda enviar los resultados del cálculo a otro programa de postproceso diferente.

Como regla habitual, no será necesario modificar este tercer grupo de unidades por parte del usuario³ pero se permite acceso a ellas por las razones indicadas anteriormente.

(¹) Este criterio tiene dos aspectos, el primero sería de orden psicológico y de costumbre en el usuario, y el segundo iría en relación al calibrado correcto de las tolerancias en la modelación geométrica y en la generación de malla.

(²) Cada uno de estos grupos se expresaría mediante una ventana distinta en el programa.

(³) Estas consideraciones se explicarán claramente en la ventana de unidades del programa

2.6.4 Unidades en el postproceso

El hecho de que la mayoría de programas de análisis no incluyan el concepto de unidades en el conjunto de sus datos¹, significa que éstas no podrán venir dadas al postproceso desde el módulo de análisis mismo. El programa de análisis normalmente no conoce en que unidades están expresados los datos de entrada al cálculo. Así, le será imposible transmitirlas al posproceso. Será necesario definir algún sistema para hacer esta transmisión de unidades desde el pre hacia el postproceso. El origen de la problemática viene determinado por varios factores:

• El postproceso lee los resultados como simples magnitudes escalares, vectoriales o tensoriales. Desconoce el significado físico de cada uno de estos resultados y, por tanto, no puede deducir automáticamente de que magnitud física tratan.
• La libertad que se da al usuario de ejecutar el análisis tanto desde el interior del sistema como mediante la escritura de un fichero, supone la pérdida de control desde el sistema respecto a la vinculación de un postproceso con un preproceso.

El método propuesto para solventar esta problemática es el siguiente: cuando se define el tipo de problema, se deben explicitar las magnitudes físicas para cada uno de los resultados. La información podría introducirse en el mismo fichero donde se han definido las unidades del tipo de problema en el preproceso (ver página actual). La información se describiría según la Tabla 2.

Tabla. 2 Magnitudes de unidades para los resultados en postproceso

Resultados		Magnitudes
RESULT:	desplazamientos	length
RESULT:	tensiones nodales	pressure
…

El preproceso podría leer estos datos, cambiar las unidades según las definidas en el grupo tercero y escribir un fichero adicional, paralelo al de cálculo, donde se escribiría esta información. Una vez realizado el cálculo, el postproceso leería este fichero al mismo tiempo que los resultados. La visualización de resultados podría hacerse en las unidades del grupo de preferencias segundo, escalando los valores de resultados según sea necesario y visualizando en la pantalla gráfica las unidades en que se representa cada resultado.

La razón de que deba escribirse este fichero adicional de comunicación entre el pre y el postproceso, radica en el hecho de que deben almacenarse las unidades que efectivamente se usaron en el momento de escribir el fichero de datos para el análisis y no las que haya en un momento posterior en el programa.

(¹) Hay un número elevado de programas de análisis que únicamente requieren que los datos de entrada estén expresados en unidades coherentes entre ellas, sin importar cual es su nombre y por tanto, sin poder realizar transformaciones internas de unas unidades a otras.

2.7 APORTACIONES A LA ORGANIZACIÓN DE GiD

Se han definido un conjunto de criterios para el diseño del sistema y para una interacción correcta con el usuario. Estos criterios están destinados tanto a minimizar la curva de aprendizaje del usuario como para permitirle un uso cómodo y eficiente de toda la funcionalidad que se ofrece.

Se ha subdividido el conjunto de información que se requiere para una análisis genérico en varios subgrupos lógicos y funcionales que más tarde, en el Capítulo 5, permitirán introducir el concepto de adaptación.

Se ha descrito con detalle la manera de tratar los sistemas de unidades dentro de GiD y su relación con los datos que se transmitarán a los códigos de análisis.

2.8 CONCLUSIONES A LA ORGANIZACIÓN DE GiD

El diseño de un sistema que interactúe con el usuario y permita realizarle un conjunto de operaciones complejas, puede considerarse una tarea difícil. En concreto, ha sido necesario definir un conjunto de criterios que rijan el comportamiento general del sistema y que lo acerquen a ese óptimo de facilidad y eficiencia de uso.

Asimismo, la subdivisión de la información interna en subgrupos y la forma de uso de las unidades permite la implementación práctica de un sistema que se pueda adaptar a cualquier tipo de análisis y posibilita la introducción del concepto de adaptación en capítulos posteriores.

3 Tratamiento y modelización geométrica

En este capítulo se describirán las estructuras de tratamiento y almacenamiento de la información geométrica así como su visualización, creación y manipulación por parte del usuario. También se describirán algunos de los algoritmos usados para la manipulación de las entidades geométricas.

Debe entenderse que parte de la funcionalidad de las técnicas y algoritmos que van a ser descritos, está compartida con cualquier programa genérico que trate con geometrías, como puede ser un sistema de CAD o CAM. En el otro sentido, hay un conjunto amplio de algoritmos que están especializados en preparar a las entidades geométricas para su mallado posterior. Normalmente, se usa el concepto CAE, para referirse al uso de programas de tratamiento geométrico para realizar análisis posteriores por ordenador.

3.1 DEFINICIONES

Para introducir toda la temática de la modelización geométrica, será interesante ofrecer una serie de definiciones sobre los conceptos básicos.

Entidad geométrica: Figura geométrica que puede ser tanto un punto como una línea, superficie o volumen, que tiene una definición matemática y una representación gráfica. Estas entidades pueden llegar a tener una malla asociada o ser solamente auxiliares.

Geometría booleana de sólidos: Tratamiento geométrico basado en sólidos, o sea un conjunto de entidades geométricas que definen un volumen en el espacio. Típicamente se incluye en esta geometría una serie de operaciones que afectan a uno o varios de estos sólidos como: uniones, intersecciones, extrusiones, etc.

Malla: Conjunto de elementos geométricos simples como triángulos o tetraedros, unidos entre ellos de manera que comparten sus nodos en el caso de las mallas conformes y que son independientes en el caso de las no conformes, que sirven para describir de manera sencilla la geometría de un continuo. Las mallas se usan para describir el modelo en la mayor parte de los métodos numéricos de análisis por ordenador.

3.2 DESCRIPCIÓN Y COMPARACIÓN DE DIVERSOS MODELOS GEOMÉTRICOS

Para describir la jerarquía de entidades geométricas y como están tratadas y vinculadas, es conveniente hacer una comparación con los sistemas que se usan habitualmente en programas existentes. Uno de los factores que se tendrán más en cuenta en la comparación será la generación posterior de la malla, que es el objetivo principal del tratamiento de la geometría en nuestro caso.

3.2.1 CAD tradicional

En un sistema de CAD tradicional, cada entidad es independiente de las entidades contiguas. Esto significa que dos líneas unidas por uno de sus extremos, en realidad sólo tienen iguales las coordenadas de posición de esos puntos extremos. No hay ningún tipo de vinculación topológica entre ellas. Por tanto se podrían describir como entidades con información geométrica pero no topológica¹.

Figura 3: Si no hay información topológica de asociación de entidades, no se puede distinguir entre estas dos configuraciones.

Para poder generar una malla sobre las entidades geométricas, este modelo presenta serios inconvenientes. La malla debe tener relaciones topológicas claras entre un elemento y sus vecinos. En función de si comparten o no comparten los mismos nodos dos elementos vecinos, se estarán modelando dos comportamientos distintos. Ello implica que en el momento de generar la malla, debe realizarse algún tipo de operación de relación de entidades entre ellas. Este tipo de operaciones pueden llegar a ser muy costosas y difíciles.

En la Figura 3, se pueden ver dos configuraciones que deberían generar dos mallas completamente distintas. Para ello hay que mantener información de estas conectividades.

(¹) Es habitual en los programas de CAD, que tengan una serie de utilidades de creación de entidades a partir de datos geométricos de otras, como podría ser un punto final de línea o una intersección. Ello no significa que se mantenga esta información de relación en la base de datos del sistema.

3.2.2 CAD con geometría booleana de sólidos

En otros sistemas de CAD el tratamiento de la geometría se hace mediante la geometría booleana de sólidos. En este caso sí hay información topológica que relaciona todas las entidades geométricas que conforman un sólido o volumen. Esta información es imprescindible para poder realizar las operaciones de sólido sobre todo el conjunto de entidades. En concreto se tendrá información sobre todas las superficies que conforman un sólido y sobre las relaciones de vecindad entre ellas.

Figura 4: Una operación típica en la geometría booleana es la unión entre dos sólidos formando un tercer sólido.

La información de sólido sí que puede ser interesante para una generación de malla posterior. En concreto, si los modelos a tratar son únicamente volúmenes, se tendría toda la información topológica necesaria. El problema aparece cuando en el modelo existen también superficies o líneas aisladas que podrían ser usadas para simular elementos de lámina o de barra. También ocurre que en modelos de análisis numérico puede ser interesante definir elementos especiales como elementos de junta, de contacto entre sólidos, etc. En estos casos es necesario extender la información topológica de solidos al conjunto de todas las entidades.

3.2.3 Superficies como contenedoras de superficies

Otra alternativa para almacenar información de relación es la de usar unos tipos especiales de entidades que actúan como contenedores de un conjunto de superficies. El usuario o el mismo programa definiría estas polisuperficies o supersuperficies a partir de un cierto conjunto de superficies y así se crearía información topológica.

Esta información también puede usarse para una generación de malla pero se mantiene la problemática de que no cubre toda la casuística necesaria.

3.2.4 Modelo jerárquico de entidades geométricas

En este modelo, se define una jerarquía de entidades donde hay cuatro diferentes niveles:

• Primer nivel: puntos.

• Segundo nivel: líneas.

• Tercer nivel: superficies.

• Cuarto nivel: volúmenes.

Toda entidad del nivel segundo al cuarto estará apoyada en sus extremos en entidades de nivel inferior. A continuación se describen estas relaciones:

Figura 5: En la figura se pueden apreciar las relaciones jerárquicas entre diferentes entidades geométricas.

• Los dos extremos de una línea contendrán una entidad punto cada uno. Este punto será también una entidad geométrica y será el mismo punto para dos líneas que compartan un extremo.

• El contorno de una superficie estará formado por un conjunto cerrado de entidades línea. En el caso de que la superficie tenga agujeros, serán varios conjuntos cerrados. Dos superficies que estén juntas por un lado, compartirán la misma línea o líneas.

• Para el caso de volúmenes, se opera de manera similar a la de las superficies.

En la Figuras 5 y 6, pueden observarse algunos ejemplos de estas relaciones.

De esta manera se consigue tener todas las entidades relacionadas entre ellas. O sea, que el mismo usuario indica qué es lo que debe estar unido y lo que no. De alguna manera, al mismo tiempo de estar definiendo geométricamente el modelo, lo está definiendo topológicamente.

Como inconveniente a este modelo, se puede citar la mayor dificultad en las operaciones sobre la geometría tanto a nivel de usuario como a nivel del desarrollo del código y de los algoritmos. A nivel de usuario, este tiene que tener en cuenta en todo momento la forma con la que están unidas las diferentes entidades y si algunas entidades deben colapsarse con otras al realizar algunas operaciones. A nivel del desarrollo del código, aparecen una serie de complicaciones debidas al hecho de que a las operaciones geométricas hay que añadir un control de la información topológica. En algunos casos esta dificultad añadida requiere complicados algoritmos de búsqueda y criterios de comparación.

La gran ventaja de este modelo es que la geometría reproduce de manera natural la topología de la malla. Por tanto, una buena definición del modelo geométrico, deja a la malla perfectamente definida a nivel de conexión de elementos y permite mezclar con facilidad elementos de diferente tipo como pueden ser de barra, de lámina, de contacto o de sólido.

Figura 6: Jerarquía de entidades para un volumen.

3.3 MODELO ESCOGIDO EN LA TESIS

A partir de la comparación del apartado anterior, se escogió el modelo jerárquico por ser el que mejor se adapta a las necesidades de un preprocesador para análisis por métodos numéricos. Se debe tener en cuenta que la obtención de una malla para el análisis es el objetivo principal buscado. Las dificultades adicionales provenientes de la manipulación geométrica pueden sortearse mediante la importación de geometrías desde un CAD tradicional.

Por otra parte, unos modelos no son excluyentes de otros. Por tanto podría pensarse en adaptar a este modelo la geometría booleana de sólidos, o el concepto de polisuperficies para permitir un tratamiento más potente de la geometría sin perder el control necesario para la generación de la malla.

Para definir el modelo usado para tratar la geometría, tan importante como describir las relaciones topológicas entre las entidades, es exponer que tipos de entidades geométricas se van a usar y que descripción matemática tienen.

Para su exposición, se van a dividir en cuatro grupos dependientes de su dimensión interna que son: puntos, líneas, superficies y volúmenes.

3.4 ENTIDADES PUNTO

Una entidad punto se define geométricamente por sus coordenadas. Independientemente de los múltiples medios que existen para crear un punto, como pueden ser ejes locales, coordenadas relativas, sistemas de coordenadas cilíndricos o esféricos, etc. el almacenamiento interno se realiza mediante sus coordenadas en el sistema global de ejes coordenados.

Figura 7: Un punto en la geometría fuerza un nodo en la malla.

El punto servirá como base para las subsiguientes líneas y para forzar que se genere un nodo en la malla en el posición geométrica de dicho punto. También tendrá interés para asignar condiciones de contorno sobre puntos. En la Figura 7, puede verse un ejemplo de este caso.

3.5 ENTIDADES LÍNEA

Las entidades línea se definen topológicamente como la unión de dos entidades punto. En caso de ser una línea cerrada, el punto inicial y el final coincidirán.

Los diferentes tipos de entidades geométricas, que se diferencian por su descripción matemática y que se se incorporan en este trabajo son:

• Líneas rectas Aunque este tipo de entidad puede representarse como una línea NURBS lineal, interesa considerarla como una entidad geométrica diferenciada. La ventaja de hacer un caso especial con ella es que se almacena menos información con la representación simple y se simplifican algunos algoritmos.
• Arcos De igual manera que en las líneas rectas, éstos pueden representarse con una NURBS racional pero interesa más mantener la representación clásica de un arco bidimensional más una matriz de transformación asociada. La ventaja de este segundo esquema consiste en que se consigue una parametrización con el parámetro arco, lo cual tiene indudables ventajas en la generación de malla.
• Polilíneas Se entiende como polilínea a un conjunto ordenado de varias líneas simples. Su mayor ventaja reside en que simplifica el tratamiento de un conjunto grande de pequeñas líneas una detrás de otra a la hora de seleccionarlas o aplicarles cualquier operación.
• Líneas NURBS Es la entidad más potente y más versátil de todas las descritas. Con ella se pueden describir exactamente polígonos de cualquier grado. También permite describir cualquier tipo de línea racional como círculos, arcos, elipses, etc. Pero lo que la hace sin duda más potente es el hecho de que permite crear curvas a trozos, o sea, que la curva final es la unión de un conjunto de polinomios o curvas racionales con suficiente continuidad en los puntos de unión. Todas estas propiedades han convertido a las líneas NURBS en el estándar de uso en cualquier sistema de CADG¹ tridimensional. La descripción detallada de la matemática asociada a este tipo de líneas puede encontrarse en el apéndice 8.

Ejemplos de estas líneas pueden verse en la Figura 8.

Las líneas en un programa tridimensional pueden entenderse de varias maneras. En concreto su uso puede ser:

• La entidad definidora de elementos de barra.
• El contorno y medio de unión o no unión de las superficies.
• La entidad donde se asignan condiciones de contorno que más tarde se transferirán a sus elementos o nodos.
• Zona de asignación de tamaños de malla para realizar concentraciones de elementos.

Figura 8: Diferentes tipos de líneas implementadas.

(¹) El acrónimo CADG proviene de Computer Aided Design Geometric, que se traduce por diseño geométrico asistido por ordenador.

3.6 ENTIDADES SUPERFICIE

Las superficies pueden definirse de manera equivalente a las líneas en cuanto a criterios topológicos. En concreto, una superficie tiene como contorno a un conjunto de líneas y puede ser ella misma contorno de una entidad volumen.

Figura 9: Diversos tipos de superficies implementadas.

En cuanto a su descripción matemática, podemos diferenciar los siguientes tipos de superficies implementados en la tesis:

• Superficie plana únicamente definida por las líneas de su contorno y por la descripción del plano. Se simplifican algunas operaciones como la generación de malla de superficie. Su inconveniente es que, al no ser paramétrica, no se tienen direcciones preferentes para operaciones tales como generación de malla estructurada o cálculo automático de ejes locales.
• Superficie Coon es un interpolante de las líneas de contorno. Por tanto, la única información geométrica viene dada por esas líneas. Esta propiedad simplifica algunas operaciones como las distorsiones de geometría ya que sólo es necesario realizarlas sobre el contorno. Su mayor inconveniente reside en que no se puede modelizar formas complejas con ellas.
• Superficie NURBS De la misma manera que sus homólogas en líneas, es la entidad más potente para realizar modelación geométrica y por tanto, se han convertido en el estándar de todo tipo de aplicaciones tridimensionales. Permiten modelizar cónicas así como superficies polinómicas a trozos. Otra de sus propiedades reside en que es posible modificar su forma actuando sobre su polígono de control. Con ello se puede realizar una modificación interactiva de forma en la fase de diseño. Una descripción detallada de este tipo de superficies puede encontrarse en el apéndice 8.
• Superficie NURBS recortada¹ Es una extensión de la anterior en la que se permite que las líneas de contorno no se correspondan con el contorno natural de la NURBS. El único requerimiento es que pertenezcan al interior de la superficie. Típicamente se obtienen a partir de intersecciones de superficies simples.
• Superficie-malla es un tipo de superficie cuya información geométrica viene dada por una malla de triángulos o cuadriláteros. Este tipo de entidades se describirá en el Apartado 3.10.7.

Ejemplos de estas superficies pueden verse en la Figura 9.

Su uso puede entenderse de igual manera que el de las líneas sin más que cambiar elementos de barra por elementos de superficie en el Apartado 3.5.

(¹) Es habitual referirse a ellas como superficies NURBS trimadas.

3.7 VOLÚMENES

Figura 10: En la figura se puede apreciar un volumen cuyo contorno exterior es un cono.

Se define volumen como el espacio comprendido en el interior de un conjunto cerrado de superficies que definen su contorno. Por extensión, el volumen puede tener agujeros que son conjuntos cerrados de superficies que definen un espacio vacío en su interior.

La única información geométrica que los define son el conjunto de superficies que actúan como su contorno. En esta entidad no es necesaria información geométrica interior. Un ejemplo de la representación de un volumen puede verse en la Figura 10.

3.8 CONTACTOS

Hay algunos tipos de simulaciones numéricas que requieren elementos especiales de unión entre diferentes partes del modelo. Estos contactos pueden ser de elementos degenerados o pueden ser de movimiento de sólido rígido. En ambos casos debe establecerse una relación nodo a nodo en la malla que se va a generar.

En los contactos de elementos degenerados, típicamente se duplica una línea en 2D o una superficie en 3D. Ambas entidades ocupan la misma posición geométrica y pueden ser la separación de dos volúmenes que están uno al lado del otro o pueden ser el contacto entre un volumen y el exterior (en 2D, lo mismo para superficies). El contacto que se crea entre ambas entidades produce, al mallar, elementos degenerados que pueden ser tanto de línea, unión nodo a nodo, como cuadriláteros, hexaedros o prismas de base triangular. En todos los casos, una de sus dimensiones será 0.

Figura 11: Ejemplo de un contacto entre volúmenes

Los contactos con movimiento de sólido rígido, se pueden crear entre dos entidades con dimensiones geométricas equivalentes, pero con un movimiento de sólido rígido entre ellas. En este caso, también hay una relación nodo a nodo entre ambas superficies y se pueden crear elementos de contacto equivalentes al caso anterior. Típicamente se crean para forzar que se genere exactamente la misma malla en dos zonas del modelo.

Se han implementado dentro del sistema todas estas posibilidades de contacto mediante la inclusión de una entidad considerada como de superficie, para problemas 2D y como de volumen, para problemas 3D, que permite forzar esta coincidencia de nodos y la creación de los elementos descritos.

Entre la multitud de casos donde puede ser interesante usar contactos se pueden citar:

• Contactos entre dos piezas mecánicas para evitar que una penetre en el interior de la otra.
• Modelización de fenómenos de rozamiento.
• Modelización de juntas en mecánica del suelo.
• Transmisión térmica entre sólidos o entre un sólido y el exterior.
• Simetría de rotación en análisis de fluidos.

La representación gráfica de un contacto puede verse en la Figura 11.

3.9 IMPORTACIÓN DE MODELOS GEOMÉTRICOS

Una capacidad absolutamente imprescindible en cualquier programa de modelado sólido o de generación de malla, es la de poder intercambiar modelos entre diferentes programas. Esta necesidad aparece en los programas de modelado geométrico para permitir el intercambio de información entre personas que no usen el mismo programa. En el caso de un preprocesador para análisis numérico, esta necesidad se hace aún más imprescindible debido al hecho de que, típicamente, el modelo que se quiere analizar ya ha sido creado para otras utilidades como el diseño, la mecanización o la definición de dimensiones.

A la dificultad inherente a todo tipo de importación geométrica se suma, por tanto, la dificultad añadida de que el modelo en muchos casos no se ha diseñado pensando en el análisis. Ello implica que la geometría puede ser no apta, a priori, para realizar una generación de malla sobre ella.

Otras dificultades provienen del hecho de que las entidades geométricas que se definen el el CAD no han sido pensadas para generar malla sobre ellas. Además la escritura, transformación y lectura es un proceso que introduce pequeños¹ errores en la información geométrica que es necesario corregir.

Consecuencia de todo esto es que resulta conveniente dividir el proceso de importación en dos fases diferenciadas:

• La importación propiamente dicha. Consiste en la transformación de las entidades del fichero de intercambio al formato de la base interna de datos.
• La corrección y adaptación de la geometría para convertirla en apta para la generación de malla sobre ella.

En los siguientes apartados se van a describir ambas fases. La primera de una manera genérica y la segunda con mayor precisión, pues es la que tiene una relación más alta con el concepto de preproceso.

Figura 12: Diagrama de flujo del proceso de importación de modelos geométricos.

(¹) No es nada extraño encontrar discrepancias en la geometría importada que provoquen errores en la tercera e incluso en la segunda cifra significativa.

3.9.1 Importación y ficheros de intercambio

Cada programa de modelado geométrico tiene su propio formato de base da datos interna. Esto se hace así para facilitar y optimizar el almacenado de las entidades que maneja el código. El inconveniente de este sistema de organización de datos es que unos programas no pueden entender la base de datos de los otros. Ello conlleva a definir unos formatos de fichero de intercambio, aceptados por todos los programas, que son los que se van a usar para importar y exportar los modelos. A lo largo de los años, una serie de instituciones y entidades han ido definiendo varios de estos estándares y algunos de ellos son aceptados y usados por la mayoría de CAD y programas de tratamiento geométrico. En este apartado se van a describir brevemente tres de ellos, estudiados e implementados en el programa desarrollado en la tesis.

1. Formato DXF Este formato se creó y se usa mayoritariamente para el intercambio de geometrías bidimensionales. Soporta la definición de entidades tridimensionales pero de una manera un tanto precaria. En concreto, no soporta las superficies NURBS.
2. Formato IGES Este formato es inherentemente tridimensional y se puede considerar como el estándar para todos los programas que trabajan con geometrías complejas en el espacio. Soporta gran cantidad de entidades diferentes, lo cual es, en algunos casos, un inconveniente por la gran complejidad que entraña. Soporta las líneas y superficies NURBS.
3. Formato VDA Es un formato tridimensional simplificado en el cual todas las superficies tienen que convertirse a unas entidades básicas que tienen forma analítica. Su mayor ventaja es la simplicidad, lo cual evita muchos errores en el traspaso. Su inconveniente es que, al pasar de un sistema a otro, ya no se tiene exactamente la misma entidad que en el original. No soporta las superficies NURBS.

El código encargado de leer estos formatos debe ser capaz de transformar todas las entidades que no están directamente soportadas en el sistema a otras que si lo estén. En algunos casos, ello implica una aproximación, pero si se tiene como entidad básica a las NURBS, casi cualquier entidad soportada en esos formatos puede ser descrita exactamente como una superficie de este tipo.

3.9.2 Filtros de corrección de la geometría importada

El proceso descrito en el apartado anterior, acaba cuando se tienen dentro del sistema las mismas entidades que estaban descritas en el fichero de intercambio. Pero el proceso no acaba aquí. Es necesario convertir ese conjunto de entidades en un modelo geométrico apto para la posterior generación de malla.

Gran parte de esta mejora del modelo puede realizarse de forma automática, o sea, sin intervención del usuario. En algunos casos es necesario que se corrijan manualmente algunos de los problemas que puedan aparecer. Esto pasa especialmente cuando el modelo no se ha diseñado pensando en que se debía utilizar para una simulación numérica o CAE. Para estos casos, un buen programa de preproceso debe tener, además de los filtros de corrección, un buen número de herramientas de manipulación geométrica que permitan al usuario realizar dichas mejoras y adaptaciones sobre el modelo.

En los siguientes apartados se describirán una serie de filtros, o algoritmos de corrección y mejora, que se aplican en cadena a un modelo geométrico proveniente de un fichero de intercambio.

3.9.3 Colapsado de entidades

Los formatos de intercambio mencionados incluyen la descripción matemática de las entidades, pero no describen la vinculación topológica que hay entre ellas. Además, es habitual que haya discrepancias entre valores que deberían ser equivalentes para una correcta conexión. Para solventar este problema es necesario definir una tolerancia que, para distancias menores que ella, se considere que dos entidades deban convertirse en una. Esta tolerancia pueda ser dada por el usuario o puede calcularse automáticamente como un determinado porcentaje del tamaño total del modelo.

Figura 13: El proceso de colapsado reduce las cuatro líneas iniciales a tres.

Dada una tolerancia ${\textstyle \epsilon }$, el criterio de colapso según entidades será el siguiente:

• Dos puntos se convertirán en uno si se cumple que:

$\|{\vec {P}}_{1}-{\vec {P}}_{2}\|{<\epsilon }$

• Dos líneas se colapsarán si comparten los puntos finales y la distancia máxima entre ellas es menor que ${\textstyle \epsilon }$. También se realizará el colapso si una de ellas pertenece al interior de la otra. O sea, que sus puntos finales ${\textstyle {\vec {P}}_{i}}$ cumplan que

$\mathrm {dist} ({\vec {P}}_{i},L)<\epsilon$

además del criterio anterior. En este último caso, el producto del colapso no será una única línea sino varias.

• Dos superficies se colapsarán si comparten las mismas líneas de contorno y cumplen un criterio de distancia máxima similar al del caso anterior.
• Dos volúmenes serán considerados equivalentes si comparten las mismas superficies de contorno.

El cálculo de la distancia máxima entre líneas o entre superficies se realiza de manera aproximada según el criterio siguiente: se escoge un número determinado de puntos ${\textstyle {\vec {P}}_{i}}$ en el interior de la línea o superficie en función de sus características geométricas¹. Se calcula la distancia punto-recta o punto-superficie ${\textstyle d_{i}}$ para cada uno de ellos. Entonces se aproxima ${\textstyle d_{max}}$ como ${\textstyle d_{max}=max(d_{i})}$.

Para el cálculo de las ${\textstyle d_{i}}$, se calcula el punto ${\textstyle {\vec {P}}_{i}'}$ como el mapping² o reposicionamiento del punto origen sobre la entidad y ${\textstyle d_{i}=\|{\vec {P}}_{i}-{\vec {P}}_{i}'\|}$. El proceso de reposicionamiento se describe en la sección 3.10.3.

Otro problema que surge al colapsar entidades, es que deja de cumplirse que las entidades de contorno de líneas o superficies pertenecen al interior de la entidad base. Este propiedad es imprescindible para que funcionen la mayoría de algoritmos posteriores, entre ellos la generación de malla. La corrección a este problema se consigue mediante el algoritmo de movimiento de puntos descrito en el Apartado 3.10.2.

Figura 14: En la figura se muestra el ahorro de elementos que supone eliminar detalles innecesarios.

El colapso también es útil para hacer desaparecer detalles excesivamente pequeños de la geometría. Si algunas líneas o superficies son menores que la tolerancia, éstas desaparecen y se unen el resto de entidades que hay alrededor. A esta operación se le llama comúnmente feature reduction o eliminación de detalles. A efectos prácticos, muchos detalles que son imprescindibles para el diseño, pueden ser obviables para el análisis. Esta capacidad puede disminuir en gran manera el número de elementos de la malla posterior. En la Figura 14, se muestra un ejemplo gráfico de tal disminución de la malla.

(¹) En general, mayor número de puntos de control y mayor grado implica una complejidad más grande. Por tanto, se escogerán más puntos cuando aumenten estos valores.

(²) Se llama mapping a la transformación geométrica sobre un punto que lo desplaza sobre otra entidad geométrica como puede ser una superficie.

3.9.4 Corrección en la lista de \emphknots Corrección en la lista de knots¹

Por definición de líneas o superficies NURBS, no puede existir en la lista de knots uno de ellos que tenga multiplicidad mayor que el orden de la entidad. Pero en algunos casos, la información contenida en el fichero de intercambio incluye tal error. En estos casos es necesaria una corrección de la lista que se basa en eliminar los knots sobrantes y eliminar los puntos y los pesos asociados. Así, dada una curva de grado ${\textstyle n}$ y ${\textstyle L}$ puntos con la siguiente lista de knots:

$k_{0},\dots ,k_{i},k_{i+1},\dots ,k_{i+n+p},\dots ,k_{L+n}\qquad {\hbox{con: }}\quad k_{i}=k_{i+1}=\dots =k_{i+n+p}$ (1) ${\vec {P}}_{1},\dots ,{\vec {P}}_{L}$ (2) $\omega _{1},\dots ,\omega _{L}$ (3)

la lista (1) quedará reducida a:

$k_{0},\dots ,k_{i},k_{i+1},\dots ,k_{i+n},k_{i+n+p+1},\dots ,k_{L+n}$

(4)

y la lista de puntos de control (2) y pesos (3) a:

${\vec {P}}_{1},\dots ,{\vec {P}}_{i+1},{\vec {P}}_{i+p+2},\dots ,{\vec {P}}_{L}$ (5) $\omega _{1},\dots ,\omega _{i+1},\omega _{i+p+2},\dots ,\omega _{L}$ (6)

Debe tenerse en cuenta que la nueva lista tendrá ${\textstyle L-p}$ puntos de control y que la curva resultante será sólo una aproximación de la curva original.

(¹) La lista de knots es un conjunto de valores reales no decrecientes pertenecientes al espacio paramétrico y que se usan para definir una NURBS. Una definición más detallada de knot se da en el Apartado 8.5.

3.9.5 Reparametrizaciones

Para poder generar malla correctamente sobre una línea o superficie, sería óptimo que la entidad estuviese parametrizada según el parámetro arco. Esta propiedad normalmente no se da en las NURBS. El problema se corrige, hasta cierto punto, con las técnicas descritas en el Apartado 4.9.2. Sin embargo, es conveniente mejorar la parametrización desde el origen para obtener mejor calidad en la generación. A estos cambios en las líneas o superficies para obtener mejor avance con el parámetro se les denomina reparametrizaciones.

Figura 15: Comparación entre derivadas de la curva antes y después de la reparametrización.

Uno de los problemas aparece cuando hay graves discrepancias entre la relación de distancias entre los nodos con la relación de distancias entre los puntos de control. La detección del problema se produce mediante una comparación entre módulos de derivadas tal que:

${\frac {1}{F}}<{\frac {\left\|\left.{\frac {d{\vec {L}}}{dt}}\right|_{t=k_{i-}}\right\|}{\left\|\left.{\frac {d{\vec {L}}}{dt}}\right|_{t=k_{i+}}\right\|}}<F$

(7)

Donde los ${\textstyle k_{i}}$ en (7) son todos los knots interiores de la lista y ${\textstyle F}$ es un parámetro máximo aceptable, cuyo valor puede situarse en 4-5 para producir parametrizaciones aceptables.

La primera técnica de corrección consiste en descomponer la curva analizada en la curva resultante de la suma del conjunto de curvas de Bézier equivalentes. Para ello se realiza una inserción múltiple de knots hasta que todos ellos tienen multiplicidad igual al orden. El proceso de inserción de knots se describe en el Apartado 8.6.3.

La curva obtenida es equivalente a la anterior tanto en su forma geométrica como en su parametrización. Dada una curva de grado ${\textstyle n}$ y ${\textstyle L}$ puntos con los siguientes valores:

$k_{0},\dots ,k_{i},\dots ,k_{L+n}$ (8) ${\vec {P}}_{1},\dots ,{\vec {P}}_{L}$ Se van insertando knots en los valores de los knots existentes en (8) hasta que todos ellos tengan multiplicidad igual al orden $\underbrace {k_{0},\dots ,k_{0}} _{n+1},\dots ,\underbrace {k_{i},\dots ,k_{i}} _{n+1},\dots ,\underbrace {k_{L+n},\dots ,k_{L+n}} _{n+1}$ (9)

El nuevo número de puntos ${\textstyle L'}$, se habrá incrementado tanto como la suma de diferencias entre multiplicidades y orden de cada uno de los knots.

Si en la nueva lista de puntos, se van escogiendo de ${\textstyle n+1}$ en ${\textstyle n+1}$, se obtienen sucesivas curvas de Bézier del mismo grado y que serán continuas en los puntos de unión entre curvas. Estas nuevas curvas no dependen de la lista de knots, por tanto, su forma no variará al trasladar cada uno de los knots. Por tanto, podemos recalcular los incrementos entre grupos de knots tal que cumplan alguna separación más favorable. Una posible modificación, basada en la parametrización chord length, consiste en que, en la nueva lista:

$k_{0},\dots ,k_{n},\dots ,k_{i},\dots ,k_{i+n},\dots ,k_{L'},\dots ,k_{L'+n}$ (10) $k_{0}=\dots =k_{n}\qquad k_{i}=\dots =k_{i+n}\qquad k_{L'}=\dots =k_{L'+n}$ (11) se calcula el incremento como: $k_{i+n+1}-k_{i+n}={\frac {l_{c}}{l_{tot}}}$ (12)

Donde ${\textstyle l_{c}}$ en (12) es la longitud de la curva de Bézier correspondiente y ${\textstyle l_{tot}}$ es la longitud total.

La curva así calculada tendrá la misma forma geométrica que la curva original pero un avance distinto con el parámetro, o sea, una parametrización distinta. En la Figura 15, puede observarse la mejora típica que se produce con la aplicación de este algoritmo.

3.9.6 Conversión a cúbica similar

En algunos casos, la corrección del apartado anterior no es suficiente para asegurar una parametrización suficientemente buena. Estos casos se dan típicamente cuando alguno de los puntos de control está repetido varias veces, o el orden de magnitud de distancias entre algunos puntos disminuye en varios factores respecto a los otros.

Figura 16: Conversión de una línea NURBS a una NURBS cúbica interpolante que la aproxima.

Para estos casos se acepta obtener una curva o superficie que no sea exactamente igual que la original pero que la aproxime lo suficiente. Debe tenerse en cuenta que la importación se realiza bajo el control de una tolerancia. Por tanto, es lícito cambiar la forma de una curva siempre que se mantenga dentro de los límites de esa tolerancia.

El proceso consiste en calcular un conjunto de puntos ${\textstyle {\vec {P}}_{i}}$ pertenecientes al interior de la curva y que definan a ésta de manera suficientemente precisa. El criterio para decidir el número de puntos ${\textstyle L_{I}}$ se obtiene de correlacionar el número de puntos de control ${\textstyle L}$, con la tolerancia aceptable. A partir de estos puntos se podrá calcular una nueva curva cúbica interpolante según las técnicas descritas en el Apartado 8.6.10. En la Figura 16, se muestra un ejemplo de la conversión. En este caso, la discrepancia entre ambas líneas es alta debido a que hay un punto con continuidad ${\textstyle C^{0}}$.

3.9.7 Unión de líneas

Dado un conjunto de líneas NURBS alineadas, es posible convertirlas en una sola línea. Si se cumplen una serie de condiciones, es provechoso tener menos líneas que definan el modelo. La razón de ello es que, en el momento de la generación de malla, el programa fuerza la inclusión de como mínimo un elemento para cada una de las entidades. En segmentos extremadamente pequeños, esto puede significar un gran incremento en la complejidad de la malla resultante.

Figura 17: Diversos criterios para aceptar o no una unión de líneas.

Los criterios para realizar la unión son los siguientes:

1. Debe existir continuidad suficiente entre los segmentos. Típicamente se considera suficiente continuidad ${\textstyle C^{1}}$.
2. Los segmentos individuales no deben tener superficies aisladas sobre ellos. El conjunto de segmentos debe ser contorno de cada una de las superficies que tienen alrededor.

Ejemplos de estas situaciones se muestran gráficamente en la Figura 17.

Para realizar la unión, las curvas deben igualarse en grado. Por tanto, se elevará el grado de todos los segmentos al del segmento con grado mayor ${\textstyle n}$. Seguidamente se unirán los puntos de control y se calculará la nueva lista de knots a partir de las originales de la siguiente manera:

$k_{A,0},\dots ,k_{A,L_{A}+n}\qquad {\hbox{curva A }}$ (13) $k_{B,0},\dots ,k_{B,L_{B}+n}\qquad {\hbox{curva B }}$ (14) entonces la nueva lista será: $\alpha k_{A,0},\dots ,\alpha k_{A,L_{A}+n},\alpha k_{A,L_{A}+n}+\beta k_{B,0},\dots ,\alpha k_{A,L_{A}+n}+\beta k_{B,L_{B}+n}$ (15) $\alpha ={\frac {l_{A}}{l_{A}+l_{B}}}\qquad \beta ={\frac {l_{B}}{l_{A}+l_{B}}}$ (16)

donde ${\textstyle l_{A}}$ y ${\textstyle l_{B}}$ en (16) son las longitudes respectivas de las líneas.

3.9.8 Subdivisión de líneas

Algunos algoritmos necesarios para operar sobre las líneas, entre ellos la generación de malla, realizan la resolución de sistemas no lineales que necesitan de la derivada de la curva para realizar los avances. Si en el interior de la curva existe una discontinuidad acusada en la derivada, estos algoritmos pueden llegar a fallar. Por esa razón, resulta interesante mantener continuidad ${\textstyle C^{1}}$ en toda la línea¹. El proceso consistirá, pues, en subdividir la curva en todos los puntos en que no se verifique la condición de continuidad adecuada.

La subdivisión se realiza mediante la inserción de knots en ${\textstyle t}$ igual al punto de discontinuidad hasta que se alcance multiplicidad igual al orden. Las nuevas curvas serán:

$k_{0},\dots ,\underbrace {k_{i},\dots ,k_{i}} _{n+1},\dots ,k_{L+n}$ se obtendrá: ${\frac {k_{0}}{k_{i}}},\dots ,{\frac {k_{i}}{k_{i}}}\qquad$ ${\hbox{curva A }}$ ${\frac {k_{i+n+1}-k_{i+n+1}}{k_{L+n}-k_{i}}},\dots ,{\frac {k_{L+n}-k_{i+n+1}}{k_{L+n}-k_{i}}}\qquad$ ${\hbox{curva B }}$

Los puntos de control simplemente se repartirán según el número de knots que correspondan a cada curva.

(¹) En general, es suficiente con que el ángulo de discontinuidad sea pequeño.

3.9.9 Reorientación del contorno de las superficies

para poder realizar la generación de malla, es imprescindible que las líneas de contorno de una superficie estén bien orientadas. Se considera una línea ${\textstyle {\vec {L}}(t)}$ bien orientada, cuando aplicando el criterio dextrógiro a su vector de avance y a la normal de la superficie ${\textstyle {\vec {N}}}$ en ese punto, obtenemos que el vector:

${\vec {V}}={\vec {N}}\times {\frac {d{\vec {L}}}{dt}}$

apunta hacia el interior de la superficie.

Figura 18: Criterio de signos y orientaciones en las lineas de recorte de una superficie trimada.

En algunos casos, la información proveniente de CAD no cumple esta norma. Es necesario realizar una comprobación y cambiar el sentido cuando sea necesario.

El problema aparece por el hecho de que no es fácil distinguir cual es el interior de la superficie. El algoritmo implementado en la tesis se aprovecha del hecho de que en las superficies no trimadas las líneas tiene que pertenecer forzosamente al contorno de la NURBS. En el caso de las trimadas, es necesario buscar alguna de las líneas que este sobre él.

Dada una línea sobre el contorno exterior de la NURBS, estará bien orientada si se cumple que:

${\frac {d{\vec {L}}(t)}{dt}}\cdot {\frac {\delta {\vec {S}}(u,v)}{\delta u}}>0\qquad {\vec {L}}\in v=0$ (17) ${\frac {d{\vec {L}}(t)}{dt}}\cdot {\frac {\delta {\vec {S}}(u,v)}{\delta u}}<0\qquad {\vec {L}}\in v=1$ (18)

Donde ${\textstyle {\vec {S}}(u,v)}$ es la superficie paramétrica. Si la línea está contenida en ${\textstyle u=0}$ ó ${\textstyle u=1}$ se emplearían expresiones similares a las de (17).

Una comprobación adicional que puede hacerse consiste en calcular la normal a la superficie en el centro calculado según los criterios del Apartado 3.10.1 y calcular la normal aproximada al contorno de líneas. Esta última se aproxima como:

${\vec {N}}_{L}=\sum _{i}{\vec {L}}(t_{i})-{\vec {C}}\times {\vec {L}}(t_{i}+\Delta t_{i})-{\vec {C}}$

(19)

Donde ${\textstyle {\vec {C}}}$ es el centro de la superficie en (19). Se considerará que las líneas están bien orientadas si:

${\vec {N}}\cdot {\vec {N}}_{L}>0$

3.9.10 Colapso de ángulos pequeños

Un problema muy habitual en los modelos geométricos importados de CAD, es que algunas superficies tiene ángulos casi nulos que imposibilitan el mallado correcto sobre ellas. En estos casos interesa colapsar ese ángulo hasta convertirlo en uno más grande que permita introducir un elemento aceptable en él.

Figura 19: Colapso de un ángulo excesivamente pequeño. Se aceptará el colapso si $L$ es mayor que un valor mínimo.

El ángulo de colapso se calculará en función de la tolerancia dada ${\textstyle \epsilon }$. En general, el criterio resultante, tal como se ve en la Figura 19, se basará en que la linea individual resultante tenga un tamaño acorde con el resto de líneas de las superficies contiguas.

3.10 ALGORITMOS DE MODELACIÓN GEOMÉTRICA

En los siguiente apartados se van a describir un conjunto de algoritmos desarrollados en la tesis que tratan esencialmente de la manipulación de la geometría para obtener diversos fines. Entre ellos:

• Proporcionar al usuario con unas herramientas de creación y manipulación geométrica para permitirle crear sus modelos y adaptarlos a la generación de malla posterior.
• Algunas de estas técnicas se utilizan indirectamente en la generación de malla, como por ejemplo el reposicionamiento de puntos sobre la superficie.
• Soporte a las técnicas de visualización de la geometría.

A efectos de lectura, cada uno de los algoritmos puede considerarse como relativamente independiente de los otros.

3.10.1 Cálculo del centro de una superficie

Para una superficie NURBS no trimada, el centro se calcula trivialmente como:

${\vec {C}}={\vec {S}}(0.5,0.5)$

En las superficies trimadas, esté cálculo se complica pues uno de los criterios imprescindibles para que muchos algoritmos funcionen, es que el centro debe estar contenido en las líneas de contorno.

Figura 20: Diferencias entre los diferentes centros de las entidades.

El proceso para la obtención del centro consiste en calcular el centro de las líneas de contorno mediante la ecuación:

${\vec {C}}={\frac {\sum _{i,j}{\vec {L}}_{j}(t_{i})}{i}}$

(20)

Donde ${\textstyle i}$ en (20) representa un número suficientemente grande de particiones para cada línea.

Este punto será necesario proyectarlo hacia el interior de la superficie mediante las técnicas de reposicionamiento de puntos que se describirán más adelante. En la Figura 20 se puede ver que las diferencias entre estos puntos pueden llegar a ser muy grandes.

3.10.2 Movimiento de puntos

Cuando se desea mover un punto de la geometría, si hay entidades de orden superior sobre él, estas deben ser arrastradas y deformadas acordemente. Es muy importante tener en cuenta al diseñar el algoritmo que el movimiento de puntos, líneas y superficies debe ser coherente entre ellos para evitar que en la configuración final, algunas entidades que son contorno de otras, dejen de estar en su interior.

Figura 21: El movimiento de un punto provoca distorsiones en las entidades vecinas.

Para ello se calcula una esfera que sea contenedora de todas las entidades que están conectadas directamente al punto. Si el radio de esa esfera es ${\textstyle R}$, se define una función ${\textstyle {\vec {f}}}$ en el espacio tal que:

${\vec {f}}={\frac {R-\|{\vec {P}}_{i}-{\vec {P}}_{0}\|}{R}}\Delta {\vec {P}}_{0}$

(21)

La ecuación (21) se aplicará a todas las entidades que estén contenidas en dicha esfera. Para las líneas y superficies NURBS, los ${\textstyle {\vec {P}}_{i}}$ representan los puntos de control.

3.10.3 Reposicionamiento o mapping sobre una superficie

Este algoritmo es uno de los fundamentales en la modelación geométrica y es también esencial para la generación de malla. Tal como se ha ido viendo en los apartados anteriores, gran parte de los procedimientos acaban necesitando de un reposicionamiento de puntos sobre la superficie.

En realidad, con el mismo algoritmo se resuelven dos problemas diferentes:

1. Dado un punto perteneciente a la superficie, encontrar los parámetros ${\textstyle u,v}$ que lo definen.
2. Dado un punto cercano a la superficie, encontrar el punto que sea proyección de éste sobre dicha superficie. También puede entenderse como buscar el punto

${\vec {P}}_{I}\in {\vec {S}}(u,v)\quad {\hbox{tal que }}\quad \|{\vec {P}}-{\vec {P}}_{I}\|\quad {\hbox{sea mínimo }}$

Ambos objetivos se reducen a un problema de minimización de distancias que conduce a la resolución de un sistema no lineal.

El primer caso nos da dos facilidades adicionales. Una de ellas es que se dispone del punto que es la solución a nuestro problema y por tanto es posible realizar comprobaciones de calidad de la solución. La segunda ventaja es que es más fácil llegar al punto objetivo desde un punto inicial que ya está en el interior de la superficie.

Dado un punto ${\textstyle {\vec {P}}}$ cercano a la superficie, el sistema a resolver es:

${\hbox{minimizar }}\qquad f(u,v)=\|{\vec {P}}-{\vec {S}}(u,v)\|$ (22) que es equivalente a minimizar: $f(u,v)=\|{\vec {P}}-{\vec {S}}(u,v)\|^{2}$ (23) si igualamos las derivadas parciales a 0 obtendremos: $\left[{\begin{array}{l}{\frac {\delta f(u,v)}{\delta u}}\\{\frac {\delta f(u,v)}{\delta v}}\end{array}}\right]=0$ desarrollando, ${\underset {\sim }{A}}=\left[{\begin{array}{l}2\cdot \|({\vec {P}}-{\vec {S}}(u,v)){\frac {\delta {\vec {S}}(u,v)}{\delta u}}\|^{2}\\2\cdot \|({\vec {P}}-{\vec {S}}(u,v)){\frac {\delta {\vec {S}}(u,v)}{\delta v}}\|^{2}\\\end{array}}\right]=0$ (27)

Para resolver este sistema no lineal es necesario, por una parte, tener unos buenos valores iniciales. Por otro lado, la resolución se puede ejecutar aplicando el método de Newton-Rapson con line search¹. Suponiendo unos ${\textstyle u_{0},v_{0}}$ suficientemente buenos, la resolución se realiza de la siguiente manera:

${\vec {x}}=\left[{\begin{array}{l}u_{0}\\v_{0}\end{array}}\right]$ (29) se calcula el residuo R ${\vec {R}}=-{\underset {\sim }{A}}\cdot {\vec {x}}$ (30) Se comienza a iterar ${\hbox{si }}\|{\vec {R}}\|{<\epsilon }{\hbox{ se acaba. }}$ (31) Se resuelve el sistema: ${\underset {\sim }{J}}\cdot \left[\Delta x\right]=\left[R\right]$ (32) Si no se usa line search ${\vec {x}}={\vec {x}}+\Delta {\vec {x}}$ (33) ${\vec {R}}=-{\underset {\sim }{A}}\cdot {\vec {x}}$ (34) $\|{\vec {R}}\|{<\epsilon }\qquad {\hbox{se acaba }}$ (35)

Se itera volviendo a (31) o se para cuando se alcanza un número máximo de iteraciones. Debe tenerse en cuenta que hay que controlar que en todo momento ${\textstyle u,v\in [0,1]}$.

La matriz ${\textstyle {\underset {\sim }{J}}}$ en (32), es el jacobiano de la transformación y su expresión es:

${\underset {\sim }{J}}=2\left[{\begin{array}{cc}\left\|\left({\frac {\delta {\vec {S}}(u,v)}{\delta u}}\right)^{2}({\vec {S}}(u,v)-{\vec {P}}){\frac {\delta ^{2}{\vec {S}}(u,v)}{\delta u^{2}}}\right\|^{2}&\left\|{\frac {\delta {\vec {S}}(u,v)}{\delta u}}{\frac {\delta {\vec {S}}(u,v)}{\delta v}}({\vec {S}}(u,v)-{\vec {P}}){\frac {\delta ^{2}{\vec {S}}(u,v)}{\delta u\delta v}}\right\|^{2}\\\left\|{\frac {\delta {\vec {S}}(u,v)}{\delta u}}{\frac {\delta {\vec {S}}(u,v)}{\delta v}}({\vec {S}}(u,v)-{\vec {P}}){\frac {\delta ^{2}{\vec {S}}(u,v)}{\delta u\delta v}}\right\|^{2}&\left\|\left({\frac {\delta {\vec {S}}(u,v)}{\delta v}}\right)^{2}({\vec {S}}(u,v)-{\vec {P}}){\frac {\delta ^{2}{\vec {S}}(u,v)}{\delta v^{2}}}\right\|^{2}\end{array}}\right]$

(36)

El line search es una técnica para acelerar la convergencia del método mediante la disminución de iteraciones. Consiste en que al finalizar cada iteración, se acepta la dirección de avance del incremento y se realiza una aproximación lineal a la solución. Dados ${\textstyle {\underset {\sim }{A}}}$, ${\textstyle {\vec {x}}}$, ${\textstyle {\vec {R}}}$ y el recién calculado ${\textstyle \Delta {\vec {x}}}$, para obtener un ${\textstyle {\vec {x}}}$ mejorado se procede de la siguiente manera:

$\nu =1\qquad \phi _{0}=\Delta {\vec {x}}\cdot {\vec {R}}\qquad {\vec {x}}_{0}={\vec {x}}$ (37) ${\vec {x}}={\vec {x}}_{0}+\nu \cdot \Delta {\vec {x}}$ (38) ${\vec {R}}=-{\underset {\sim }{A}}\cdot {\vec {x}}$ (39) $\phi =\Delta {\vec {x}}\cdot {\vec {R}}$ (40)

$\left\|{\frac {\phi }{\phi _{0}}}\right\|<\epsilon \qquad {\hbox{se acaba. }}$ (41) se actualiza $\nu$ $\nu =\nu {\frac {\phi _{0}}{\phi _{0}-\phi }}$

(42)

Y se itera volviendo a (38). Típicamente se realizan 3 ó 4 iteraciones que mejoran el acercamiento a la solución.

Los valores iniciales pueden venir del proceso de cálculo que llama a esta función. Un ejemplo de ello sería que cuando se realiza la generación de malla por el métode de avance frontal en el espacio y se debe reposicionar el punto, se pueden aproximar los valores de ${\textstyle u}$ y ${\textstyle v}$ a partir de los ${\textstyle u_{i},v_{i}}$ de los puntos del frente usados.

Generalmente es imprescindible realizar un acercamiento a la solución mediante algún método simple. En el caso de superficies, el método usado, dada una granularidad ${\textstyle N}$ que puede ser del orden de 7, consiste en:

$u_{0}=v_{0}=0\qquad u_{1}=v_{1}=1$ (43) ${\vec {P}}_{i,j}={\vec {S}}(u=u_{0}+{\frac {(i-1)}{N}}(u_{1}-u_{0}),v=v_{0}+{\frac {(i-1)}{N}}(v_{1}-v_{0}))\qquad i,j=1,\dots ,N+1$ (44) ${\vec {P}}_{n}={\vec {S}}(u_{n},v_{n})\quad {\hbox{tal que }}\quad \|{\vec {P}}_{n}-{\vec {P}}_{i,j}\|\quad {\hbox{sea mínima }}$ (45) $u_{0}=u_{n}-{\frac {1}{N}}\qquad u_{1}=u_{n}+{\frac {1}{N}}$ (46) $v_{0}=v_{n}-{\frac {1}{N}}\qquad v_{1}=v_{n}+{\frac {1}{N}}$ (47) ${\hbox{Se itera volviendo a (44). }}$

Este proceso se realiza varias veces y puede combinarse con la resolución por Newton-Rapson pasando repetidamente de una metodología a la otra.

(¹) Búsqueda según una determinada dirección de avance.

3.10.4 Intersecciones de líneas y líneas-superficies

La intersección de líneas en el espacio es, en realidad, la búsqueda de ${\textstyle {\vec {P}}_{1}\in {\vec {L}}_{1}(t)}$ y ${\textstyle {\vec {P}}_{2}\in {\vec {L}}_{2}(t)}$, tal que ${\textstyle \|{\vec {P}}_{1}-{\vec {P}}_{2}\|}$ sea mínima. La solución a este problema consiste en resolver un sistema de ecuaciones no lineal similar al del Apartado 3.10.3. La particularidad que debe tenerse en cuenta en este caso, es que dos líneas pueden intersectarse en varios puntos y algunos de ellos puede ser ya un punto final. Para estos casos hay que partir de una cierta distancia desde los puntos ya conectados cuando se aplica el algoritmo de aproximación de la solución.

Figura 22: Zona válida para realizar la aproximación en la búsqueda de la intersección de dos líneas en el caso de que ya estén unidas.

En la Figura 22, se dibuja la zona válida para aplicar el primer algoritmo básico de aproximación.

La intersección entre líneas y superficies se resuelve de manera similar. Cuando la línea y la superficie no intersectan, se puede aplicar la misma técnica para extender la línea hasta que entre en contacto con la superficie.

3.10.5 Intersección de superficies

Figura 23: En la figura (A), la búsqueda de la línea de intersección se realiza mediante la contraposición de cuadrículas de puntos. En la (B), se parte de una intersección línea de contorno-superficie.

La intersección de superficies consiste en calcular todas las líneas que pertenecen a la unión entre ambas entidades. Estas líneas serán una aproximación numérica a la intersección real. El proceso se puede subdividir en dos partes.

1. Encontrar algún punto de partida para todas las líneas que forman parte de dicha intersección.
2. Calcular la expresión matemática de estas líneas como aproximación a las reales.

Para realizar el primer paso, primero se deben encontrar las intersecciones de todas las líneas de contorno de cada una de las superficies con la otra mediante las técnicas descritas previamente. Estos puntos serán origen de algunas de estas líneas. El resto de las líneas deben encontrarse con algún algoritmo de búsqueda. El algoritmo implementado se basa en la tecnología de Octtrees descrita en el Apartado 4.12.1. Se define una cuadrícula de puntos suficientemente densa para cada una de las superficies y se van comparando los de una con los de la otra mediante la búsqueda rápida en el árbol de puntos. Cuando se encuentran puntos suficientemente cercanos, se considera la zona como de posible paso de una línea de intersección. Para las superficies trimadas, deben excluirse las zonas que pertenecen al exterior de las líneas de recorte.

Conociendo algún punto perteneciente a cada una de las líneas, el cálculo efectivo de éstas se realiza mediante un algoritmo iterativo. Dado ${\textstyle {\vec {P}}_{0}}$ tal que ${\textstyle {\vec {P}}_{0}\in {\vec {S}}_{1}(u,v)}$ y ${\textstyle {\vec {P}}_{0}\in {\vec {S}}_{2}(u,v)}$, se debe encontrar ${\textstyle {\vec {P}}_{1}}$ tal que:

${\vec {P}}_{1}\in S_{1}(u,v)$ (48) ${\vec {P}}_{1}\in S_{2}(u,v)$ (49) $\|{\vec {P}}_{1}-{\vec {P}}_{0}\|=\delta$ (50)

Donde ${\textstyle \delta }$ en (50) es una distancia fija que representa la precisión en la obtención de puntos para la intersección. Este planteamiento implica encontrar ${\textstyle u_{1},v_{1}}$ y ${\textstyle u_{2},v_{2}}$ tal que:

${\underset {\sim }{A}}=\left[{\begin{array}{c}S_{1,x}(u_{1},v_{1})-S_{2,x}(u_{2},v_{2})\\S_{1,y}(u_{1},v_{1})-S_{2,y}(u_{2},v_{2})\\S_{1,z}(u_{1},v_{1})-S_{2,z}(u_{2},v_{2})\\\|{\vec {S}}_{1}(u_{1},v_{1})-{\vec {P}}_{0}\|{-\delta }\end{array}}\right]=0$

(51)

La resolución de este sistema no lineal se realiza mediante las técnicas descritas en apartados anteriores. Cuando se ha llegado por ambos lados a los extremos de las líneas, que serán una línea de contorno o indicarán que la propia línea es cerrada, se detiene el proceso. Entonces se puede calcular una línea NURBS que interpole a todos estos puntos y cumpla algunos criterios de tangencia en extremos tal como se describe en el Apartado 8.6.10.

3.10.6 Visualización de superficies

Para realizar una visualización foto-realística de un modelo geométrico, existen diversas técnicas:

1. Iluminación plana. Se suministra al algoritmo de dibujo una normal por triángulo y éste pinta un color constante para cada elemento.
2. Iluminación Gouroud. En este caso se da una normal por cada vértice del triángulo. El algoritmo realiza una interpolación lineal de colores en los vértices del triángulo y luego hacia el interior.
3. Ray tracing. Para cada pixel de la pantalla, se calcula el rayo que procede del punto de vista y intersecta con ese punto de pantalla. Para cada rayo, se calcula su intersección con todos los elementos del modelo y se calculan sus propiedades lumínicas. Este modelo no está contemplado en la librería OpenGL y por tanto, no se ha usado en la tesis.

En todos los casos, además de algunas propiedades de luces y materiales, es necesario suministrar una descripción geométrica del modelo. Las librerías gráficas admiten superficies NURBS directamente pero no siempre da los resultados deseados. La solución óptima consiste en suministrar a la librería una descomposición en triángulos de la geometría original juntamente con sus normales.

Esta descomposición tiene muchos puntos en común con la generación de malla pero con la adición de algunas particularidades:

1. La conformidad de elementos entre diferentes superficies no es necesaria, ya que en la iluminación de Gouroud se promedian normales y se disimula el escalonado de los elementos.
2. Los elementos pueden tener altos niveles de distorsión. Esto permite un gran ahorro de elementos en superficies con grandes diferencias de curvatura en sus dos direcciones principales.
3. La cantidad de elementos a generar se deduce de un criterio de aproximación geométrica al modelo real. Para ello puede usarse el criterio de error cordal descrito en el Apartado 4.7.1.
4. La transición de tamaños de elemento será siempre lo más rápida posible para minimizar el número final de ellos.

En la implementación que ha dado lugar a esta tesis, la malla se genera internamente y se calculan los tamaños según un criterio de error cordal. Se permite, en los casos en que se considere apto, la generación de malla estructurada de triángulos con tamaños en ambas direcciones en función de cada una de las curvaturas principales. Finalmente, se suministran los triángulos generados juntamente con las normales nodales a la librería gráfica.

3.10.7 Las superficies-malla

Figura 24: En estas figuras se pueden apreciar dos Superficie-mallas, que provienen de una única malla de triángulos. En la figura izquierda, las superficies se representan según el modo de visualización común. En la derecha se dibuja también la malla de base.

Figura 25: La malla de esta figura es el producto de la generación sobre las entidades descritas en la Figura 24. Puede apreciarse que el tamaño de elementos es distinto en ambas.

Este tipo de entidades pueden definirse de la siguiente manera:

Superficie-malla: Es una entidad geométrica, similar a cualquier otra entidad de superficie en cuanto a su funcionalidad, cuya información geométrica viene determinada por una malla de triángulos o cuadriláteros espaciales y por unas líneas de contorno y otras interiores. Típicamente, este tipo de superficies al mallarse, proveerán una nueva malla que será muy distinta a la que les sirve de origen.

Dada una malla espacial de triángulos o cuadriláteros, el proceso de creación de una entidad de este tipo pasa por:

1. Lectura de la malla.
2. Compatibilización de normales.
3. Identificación de los contornos externos.
4. Identificación de las líneas internas, que se calculan con un criterio de ángulo máximo.
5. Subdivisión de la malla original en varias entidades de este tipo.

Una vez leída la malla, se calcula información de relaciones de entidades. En concreto, debe poder encontrarse de manera eficiente para cada uno de los elementos, la lista de sus elementos vecinos. La identificación de los contornos externos consiste en encontrar las caras de elementos que no tengan vecinos.

Los contornos internos consisten en encontrar todas las aristas entre elementos tal que el ángulo entre ambos elementos sea mayor que uno fijo. El proceso se realiza mediante avance por las aristas angulosas y subdivisión en dos líneas en caso necesario.

A partir del conjunto de aristas que forma un contorno externo o una línea interior, se interpolará una línea NURBS que cumpla criterios suficientes de tangencia. En las Figuras 24 y 25, se puede ver un ejemplo de una entidad de este tipo.

3.10.8 Creación de volúmenes

La entidad geométrica volumen se define como un conjunto cerrado de superficies que comparten sus líneas de contorno dos a dos. En el caso de que el volumen contenga agujeros, tendrá tantos subconjuntos cerrados de superficies como agujeros más el externo. Dada una selección de superficies, es necesario primero distinguir si cumplen el criterio de formar uno o varios conjuntos cerrados. Para ello, la comprobación es únicamente topológica dado que por el modelo jerárquico, las uniones entre superficies están ya definidas.

Esta comprobación se basa en que partiendo de una superficie cualquiera, se accede a las vecinas mediante el criterio de que dos superficies son vecinas si comparten la misma o mismas líneas. Cuando no se puede avanzar más, cada una de las líneas encontradas debe pertenecer exactamente a dos superficies. Se debe contemplar el caso especial de las polilíneas en las que algunas de sus sub-líneas son contorno de alguna superficie. Cuando se cierra el volumen, si quedan aún superficies en la selección, se repite el proceso para identificar los agujeros.

Mediante este mismo proceso es fácil orientar las superficies tal que todas las normales apunten hacia dentro o todas hacia fuera. Para ello, como las superficies mantienen un indicador de sentido para cada una de sus líneas, confrontando estos indicadores se puede resolver si una superficie esta en el mismo sentido o en el opuesto que la de referencia. Al acabar el proceso de búsqueda, se habrán definido los indicadores de orientación para cada una de las superficies.

3.10.8.1 Orientación de las superficies

El problema que aparece es que en este momento se tienen todas las superficies orientadas de manera coherente entre ellas, pero no se sabe si con las normales apuntando hacia dentro o hacia fuera. Para la posterior generación de la malla, es conveniente tener a las normales apuntando hacia dentro¹. El problema se resuelve calculando la magnitud con signo del volumen interior del sólido.

Por el teorema de Gauss:

$\int _{V}{\frac {\partial f}{\partial x_{i}}}dv=\int _{\partial V}f\cdot N_{i}d\Omega$

(52)

Entonces si definimos f respectivamente como ${\textstyle x_{1},x_{2}}$ y ${\textstyle x_{3}}$ se obtiene:

$\int _{V}\sum _{i=1}^{3}{\frac {\partial x_{i}}{\partial x_{i}}}dv=\int _{\partial V}\sum _{i=1}^{3}x_{i}\cdot N_{i}d\Omega$

(53)

operando:

$3\int _{V}dv=\int _{\partial V}{\vec {x}}\cdot {\vec {N}}d\Omega$

(54)

por tanto podemos obtener la magnitud orientada del volumen con la fórmula:

$V={\frac {\int _{\partial V}{\vec {x}}\cdot {\vec {N}}d\Omega }{3}}$

(55)

El signo de V nos indicará si la normal de las superficies apunta hacia dentro o hacia fuera del volumen.Para realizar esta operación se va integrando sobre cada una de las superficies y sumando sus resultados parciales.

• Para superficies planas como la normal es constante i de módulo igual al área de la superficie, se obtiene:

$\int _{\partial V}{\vec {x}}\cdot {\vec {N}}d\Omega ={\vec {x}}\cdot {\vec {N}}_{\sup }$

(56)

• Para las superficies paramétricas con parámetros u,v en el intervalo ${\textstyle \left[0,1\right]}$ como son las superficies Coon y las NURBS, se calcularía de la siguiente manera:

$\int _{\partial V}{\vec {x}}\cdot {\vec {N}}d\Omega =\int _{0}^{1}\int _{0}^{1}{\vec {x}}\cdot \left({\frac {\partial {\vec {S}}}{\partial u}}\times {\frac {\partial {\vec {S}}}{\partial v}}\right)dudv$

(57)

Donde ${\textstyle {\vec {S}}(u,v)}$ es la superficie paramétrica. La integral se resuelve numéricamente con una cuadratura de Gauss-Legendre (ver[18]) en el espacio paramétrico ${\textstyle \left[0,1\right],\left[0,1\right]}$. En caso de disponer de una tabla calculada para el espacio paramétrico ${\textstyle \left[-1,1\right],\left[-1,1\right]}$, se mantienen los mismos pesos ${\textstyle \omega _{i}}$ y se realiza una translación el los puntos de evaluación de la integral tal que:

$\xi _{i}^{N}={\frac {\xi _{i}+1}{2}}$

(58)

Para el caso de las superficies NURBS, la integración se realiza sobre las superficies de Bézier obtenidas de la descomposición de la NURBS. De esta manera, la precisión de la integral está proporcionada a la cantidad de información geométrica que contiene esta entidad geométrica.

Figura 26: Para hallar el volumen orientado de un sólido cerrado, es necesario integrar sobre las superficies.

3.10.8.2 Cálculo del centro del volumen

Para el cálculo del centro del volumen se puede usar el método elemental de asimilarlo al centroide de todos los puntos de todas las superficies.

${\vec {X}}_{c}={\frac {\sum _{i=1}^{n}{\vec {P}}_{i}}{n}}$

(59)

donde estos n puntos se obtendrían a partir de todos los puntos extremos de todas las líneas pertenecientes a su vez a todas las superficies del contorno. Este método tiene el inconveniente de que en algunos casos se concentran muchas superficies en una zona del volumen y muy pocas en otras zonas. Entonces el centro se desplaza hacia la zona en que hay más superficies y por tanto más puntos. Otro método es calcularlo mediante integración:

$\int _{V}x_{i}dv=\int _{v}{\frac {1}{2}}\partial _{i}r^{2}dv={\frac {1}{2}}\int _{\partial v}n_{i}r^{2}d\Omega$

(60)

donde ${\textstyle r^{2}=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$ y se ha aplicado ${\textstyle {gauss_{t}heorem}}$. Entonces:

$\int _{v}{\vec {x}}dv={\frac {1}{2}}\int _{\partial v}r^{2}{\vec {n}}d\Omega$

(61)

${\vec {X}}_{c}={\frac {1}{2V}}\int _{\partial v}r^{2}{\vec {n}}d\Omega$

(62)

La integración se realiza de manera similar a lo explicado en el Apartado 3.10.8.1.

(¹) En realidad, lo realmente imprescindible es saber hacia donde apuntan las normales, si hacia dentro o hacia fuera.

3.11 APORTACIONES EN LA MODELIZACIÓN GEOMÉTRICA

En este apartado se intentará distinguir entre lo que es estado del arte y conocimientos generales en contraposición a las nuevas ideas o implementaciones que se han desarrollado en la tesis. Debe entenderse como aportaciones a la resolución de problemas y creación de técnicas que no tienen una metodología comúnmente aceptada en la literatura.

3.11.1 Modelo jerárquico

Existen diversos ejemplos de sistemas que usan modelos jerárquicos. Al mismo tiempo, el concepto de información topológica asociado a la información geométrica, tiene muchos paralelos con la geometría booleana de sólidos, la cual se usa ampliamente en muchos sistemas de CAD. La aportación de esta tesis ha sido modificar el modelo jerárquico para crear un paralelo preciso entre la conexión topológica de la geometría con la conexión entre los nodos y elementos de la malla resultante. De esta manera, se pueden introducir de manera natural diversos conceptos como:

• Contactos de diversa índole.
• Creación de juntas.
• Condiciones de contorno periódicas.
• Líneas o superficies duplicadas.
• Multi-volúmenes. O sea, diversos volúmenes que ocupan posiciones superpuestas en el espacio.

3.11.2 Importación y corrección de modelos geométricos

Los modelos geométricos que provienen de CAD, han sido creados pensando en unas necesidades de diseño o de mecanización que no coinciden con las necesidades de generación de malla.

Se han desarrollado un conjunto de técnicas y algoritmos para modificar la geometría original y hacerla apta para el mallado. Entre ellos figuran los de colapsado de entidades, la corrección de los knots, las reparametrizaciones, la conversión a cúbica similar, la unión de líneas, la subdivisión de líneas, la reorientación del contorno de superficies y el colapso de ángulos pequeños. Algunas de estas correcciones solventan problemas o inexactitudes en la definición original. En otros casos, pese a que la información de base es correcta, es necesario modificarla para poder obtener en ella una malla que cumpla con los requisitos del análisis por métodos numéricos.

3.11.3 Movimiento de puntos con deformaciones

El algoritmo de movimiento de puntos arrastrando a las entidades cercanas, resuelve la problemática de mover un conjunto de entidades de manera coherente entre ellas y tal que, después del movimiento, se continúen respetando las relaciones topológicas que existían entre ellos. Así, si una entidad es contorno de otra, después del movimiento debe continuar situada geométricamente en dicho contorno.

3.11.4 Reposicionamiento de puntos sobre una superficie

Este es uno de los algoritmos fundamentales en la modelación geométrica y también se usa ampliamente en la generación de la malla. Puede entenderse, en su descripción básica, como el problema inverso de la obtención de un punto dadas sus coordenadas paramétricas, o también como un problema de minimización de distancias. La dificultad del algoritmo desarrollado estriba en que las superficies pueden presentar un muy alto grado de distorsión que dificulta o impida la resolución del sistema no lineal. Por tanto, una parte muy importante del método propuesto y descrito en el Apartado 3.10.3, está centrada en la búsqueda eficiente de unos buenos valores iniciales.

3.11.5 Intersecciones

Existe abundante literatura sobre diversas metodologías para el cálculo de intersecciones entre las diversas entidades (ver [19]). Sin embargo, éste es un problema complejo que exige buenos algoritmos aproximativos. En la presente tesis se ha optado por aprovechar todas las posibilidades que proporcionan técnicas como la de las estructuras Octtree, la interpolación mediante líneas NURBS con información de tangentes y la resolución de problemas de minimización y sistemas no lineales mediante algoritmos tipo Newton-Rapson. Esta aproximación al problema es novedosa en contraposición a las descritas en las referencias dadas [19,1,2].

3.11.6 Visualización gráfica

La técnica más ampliamente utilizada para visualizar una iluminación foto-realística de un conjunto de superficies NURBS, consiste en aprovechar las primitivas NURBS de las librerías gráficas estándares. Éstas contienen algoritmos basados en la descomposición de polígonos en triángulos y aplicados al espacio paramétrico.

En este trabajo se ha preferido aprovechar toda la tecnología implementada en generación de malla, en aproximación geométrica y en uso de superficies no conformes para cumplir este cometido. De hecho, una malla que requiera criterios de aproximación geométrica y no necesite criterios de calidad para simulación es muy similar a la malla que se necesita para visualización gráfica.

3.11.7 Superficies-malla

El uso de mallas como información geométrica base para generar nuevas mallas encima, ha sido ampliamente documentada en [20] y [11]. Sin embargo, la integración total de este tipo de entidad con el resto de entidades geométricas, dotándola de opciones similares a las otras, es una forma nueva de operar con estas superficies. Su método de mallado, totalmente integrado en el método de generación de malla en el espacio, conduce a la homogeneidad en la resolución de estos algoritmos.

3.11.8 Orientación y cálculo de volúmenes

El aprovechamiento del teorema de Gauss y de la integración numérica de superficies con cuadraturas de Gauss-Legendre para el cálculo del volumen orientado de un sólido, representa una técnica novedosa que evita la intervención humana en las decisiones relativas a la orientación de entidades.

3.12 CONCLUSIONES A LA MODELACIÓN GEOMÉTRICA

La modelación geométrica mediante líneas y superficies NURBS se ha convertido en un estándar dentro de los programas de tratamiento geométrico. Sin embargo, la aportación de la modelación jerárquica con conexión geometría-malla, aporta nuevas opciones para la simulaciones de efectos complejos tales como el tratado de contactos, juntas, volúmenes superpuestos y otros.

La parte de los desarrollos de este capítulo que introduce más novedad es, sin embargo, toda la metodología de mejora y adaptación de la geometría para la generación de malla. No se debe olvidar que el modelo geométrico no es más que un medio para la obtención de una malla que describa fielmente lo que se desea calcular. A tal fin se encaminan todas la técnicas destinadas a que el mallado sea posible, en primer lugar, y a que la malla resultante cumpla con los criterios de calidad necesarios para el análisis y simulación numérica.

4 Generación de malla

4.1 INTRODUCCIÓN

La mayor parte de los métodos numéricos actuales (elementos finitos, volúmenes finitos, diferencias finitas, integral de contorno, etc.), resuelven las ecuaciones diferenciales de gobierno de un problema sobre una malla.

Se entiende por generación de malla el proceso por el cual a partir de la geometría que define al modelo y de unas indicaciones dadas por el usuario, se procede al cálculo automático de una malla que defina adecuadamente a esa geometría con los condicionantes del programa de análisis a usar.

Para la correcta obtención de esta malla hay que tener en cuenta una serie de factores que pueden venir dados por el usuario o que pueden calcularse automáticamente:

• Tipos de elementos a generar. Unidimensionales, triángulos, cuadriláteros, tetraedros, hexaedros. De grado lineal, cuadrático, serendípitos o lagrangianos, etc.

• Selección de las entidades geométricas que deben ser malladas. Puede interesar en algunos casos obtener también la malla de contorno.

• Malla estructurada o no estructurada.

• Tamaño del elemento en cada zona.

• Transición del tamaño del elemento entre zonas.

• Comprobación de la calidad final de los elementos generados.

En siguientes apartados se comentará la problemática que origina cada uno de estos factores y las diversas soluciones que ha sido necesario desarrollar en la tesis. También se describirán los algoritmos usados en los diversos tipos de generaciones de mallas.

4.2 TIPOS DE ELEMENTOS

Según las particularidades del análisis a efectuar, puede ser necesario obtener elementos de distinto grado. Los elementos básicos, que son los que genera por defecto el código de mallado desarrollado en la tesis son:

Figura 27: Elementos de línea y triangulares

• Elementos de barra de 2 nodos (Figura 27).

• Elementos de superficie triángulos de 3 nodos o cuadriláteros de 4 (Figuras 27 y 28).

• Elementos de volumen tetraedros de 4 nodos o hexaedros de 8 (Figuras 28 y 29).

• Para contactos y elementos de junta se puede usar elementos de barra de 2 nodos, cuadriláteros de 4, hexaedros de 8 y prismas de base triangular de 6 nodos (Figura 29).

Figura 28: Elementos cuadriláteros y tetraedros

Si se desean elementos de mayor grado se pueden generar:

• Elementos de barra de 3 nodos.

• Elementos de superficie triángulos de 6 nodos o cuadriláteros de 8 ó de 9.
• Elementos de volumen tetraedros de 10 nodos.
• Elementos de volumen hexaedros de 20 ó de 27.

Figura 29: Elementos hexaédricos y prismas de base triangular.

El cálculo de los nodos adicionales para elementos de mayor grado requiere que estos nuevos nodos también estén sobre la entidad geométrica correspondiente. Por ello no es válido calcularlos simplemente como el punto medio del segmento o de la cara. Es necesario posicionarlos previamente sobre la línea o superficie a la que pertenecen o, en algunos casos, hacer un reposicionamiento¹ posterior.

Como restricciones a la elección de elementos está, por una parte, que la unión de entidades tiene que ser compatible² y por otra, que el grado de los elementos (lineal, cuadrático, o cuadrático9)³ tiene que ser el mismo para todo el modelo.

(¹) posicionamiento o mapping sobre la superficie.

(²) mallado con hexaedros no puede estar unido a uno de tetraedros. Por otra parte pueden coexistir dos volúmenes con estos tipos diferentes siempre que no haya conexión entre ellos.

(³) tipo cuadrático genera el cuadrilátero serendípito de 8 nodos y el hexaedro de 20. El tipo cuadrático9 genera el cuadrilátero lagrantiano de 9 nodos y el hexaedro de 27. Ambos generan el triángulo de 6 nodos.

4.3 ENTIDADES GEOMÉTRICAS QUE DEBEN SER MALLADAS

Dado un conjunto de entidades geométricas, el criterio por defecto es que sólo se generará malla sobre las entidades que no tengan una entidad de grado superior encima¹. Por tanto, si se tiene una superficie cuyo contorno son líneas, se obtendrá la malla de los triángulos o cuadriláteros pero no la de los elementos de barra que irían sobre las líneas. De la misma manera, si se tiene un volumen cuyo contorno son superficies, sólo se obtendrían los elementos de volumen y no los de superficie. Pero si los modelos anteriores tuvieran una línea o una superficie aislada de las anteriores, se crearía malla también sobre estas entidades.

Este comportamiento por defecto puede cambiarse para obtener, por ejemplo, además de la malla de tetraedros, la malla de triángulos que es su contorno.

(¹) Se excluye de este criterio la entidad punto por no poderse crear elementos sobre él.

4.4 MALLA ESTRUCTURADA VERSUS NO ESTRUCTURADA

Definición:Una malla es estructurada si todos sus nodos interiores tienen el mismo número de elementos que lo contienen. O sea que el número de elementos que contienen ese nodo como una de sus conectividades es constante para todos los nodos del interior del dominio.

La característica de una malla de este tipo es la regularidad. Eso significa que se puede encontrar una submalla¹ que se repite un número arbitrario de veces para formar la malla total.

Figura 30: En estas imágenes se aprecia la diferencia entre una malla no estructurada y una estructurada.

Los algoritmos de generación de malla estructurada posicionan los nodos de manera ordenada en el interior de las superficies o volúmenes. El grave inconveniente reside en que para para poder ser generada, se debe restringir la topología de las entidades geométricas que componen el modelo. De alguna manera la regularidad y la repetición de la malla resultante debe estar ya reflejada en la composición de la geometría original.

Malla no estructurada es cualquier malla que no cumpla la definición anterior. Los algoritmos de generación de malla no estructurada son más complejos y no tienen garantía absoluta de éxito pero pueden ser usados en modelos de geometría arbitrariamente compleja sin intervención del usuario en la definición de los parámetros de mallado. Es posible mezclar malla estructurada y no estructurada en un mismo modelo. En la Figura 30, se aprecian las características principales de cada tipo de malla. Más información sobre los criterios de diferenciación de tipos de malla puede encontrarse en [21].

(¹) o patrón (pattern) entendiéndose como parte de un conjunto que se va repitiendo a intervalos regulares.

4.5 MALLA CONFORME Y MALLA NO CONFORME

Se entiende por malla conforme la que todos sus elementos comparten sus nodos entre elementos vecinos y en la que cada lado de un elemento coincide con el lado de otro elemento y no con parte de ese lado. La malla no conforme es la que sus elementos son independientes entre ellos. Ejemplo de ambas mallas puede verse en la Figura 31.

En general, para modelar adecuadamente un sólido continuo, es necesario el uso de malla conforme. En algunos casos especiales, como cuando únicamente queremos discretizar un contacto geométrico, nos basta con una malla de cualquier tipo.

Figura 31: Comparación entre una malla conforme (izquierda) y no conforme (derecha).

4.6 TAMAÑO DEL ELEMENTO PARA MALLA NO ESTRUCTURADA

Definición: Se entiende por tamaño del elemento la longitud de un elemento lineal o la longitud media de una arista de un elemento de otro tipo.

En la Figura 32, se puede ver gráficamente el significado de este tamaño de elemento.

Figura 32: Tamaño del elemento para cada tipología.

Para poder generar la malla, se debe conocer el tamaño deseado de elemento en cada zona de la malla. En general, el proceso es que toda entidad geométrica puede tener información del tamaño deseado cerca de ella. Así, si se asigna un tamaño a un punto de la geometría, los elementos cercanos a ese punto tenderán a tener ese tamaño. De la misma manera, un tamaño asignado en una línea, superficie o volumen provoca que los elementos cercanos o en el interior de dicha entidad, tiendan a tener ese tamaño. Esta información puede provenir de varias fuentes:

• En el momento de generar la malla, el usuario define un tamaño medio que se asignará a todas las entidades que no tengan ningún tamaño predefinido.

• El usuario puede asignar manualmente tamaños a cualquier entidad geométrica.

• Existe una función que asigna tamaños máximos a las entidades en función de su magnitud y forma geométrica.

• Otra función asigna tamaños máximos según un criterio de curvaturas o de máximo error cordal¹ entre el modelo y la malla resultante.

• Mediante otra función se suavizan las diferencias de tamaño entre entidades cercanas para que estas diferencias no superen un determinado criterio de progresión geométrica (ver el Apartado 4.7.2). Esta función es útil después de haber asignado tamaños con algunos de los criterios anteriores.

Es interesante aplicar de manera automática el criterio de tamaños máximos por forma y el de corrección de tamaños cercanos para posibilitar de manera sencilla la consecución de una malla. En general aumentan el número final de elementos pero generan elementos de mayor corrección. En algunos casos particulares (si se necesita máximo ahorro en el número de elementos), pueden desactivarse estas funciones.

(¹) O aproximación de la máxima distancia entre la entidad geométrica y los elementos generados.

4.7 TRANSICIÓN DEL TAMAÑO DEL ELEMENTO ENTRE ZONAS

Cuando se tienen asignados diferentes tamaños a diferentes zonas de la malla, los elementos tendrán un tamaño más pequeño cerca de una de las zonas y mayor cerca de la otra. El paso de uno de los tamaños al otro puede hacerse de distintas formas. Se ha definido un parámetro adimensional de preferencia de mallado con valores entre 0 y 1 que indica si esta transición debe ser lenta o rápida. Según el tipo de análisis y el número de elementos deseado puede ser más conveniente una opción u otra. La elección de este parámetro puede influir en la calidad de la malla resultante e incluso en la posibilidad de éxito en la generación. Este parámetro influirá posteriormente en los algoritmos de generación de malla según el criterio explicado en los apartados correspondientes. La aplicación de este parámetro para una malla bidimensional de triángulos puede verse en la Figura 33.

Figura 33: El número de elementos generados cambia enormemente según el parámetro escogido.

4.7.1 Asignación de tamaños por error cordal

En algunos casos nos interesa que la malla represente la forma geométrica con una determinada precisión. Esto se consigue mediante un criterio de error cordal.

Definición: Se llama error cordal a la distancia entre la geometría original y los elementos de la malla generada. En la Figura 34, se define gráficamente dicha magnitud.

Si limitamos este error cordal a una cierta tolerancia, obtendremos una malla que aproximará tanto a la geometría como pequeña sea esta tolerancia.

Figura 34: Definición de error cordal.

Dado un error cordal ${\textstyle \epsilon _{c}}$ y una curva ${\textstyle {\vec {L}}(t)}$, el tamaño óptimo de elemento a generar sobre la curva, se calcula de la siguiente manera:

$r={\frac {\left\|{\frac {d{\vec {L}}}{dt}}\right\|^{3}}{\left\|{\frac {d{\vec {L}}}{dt}}\times {\frac {d^{2}{\vec {L}}}{dt^{2}}}\right\|}}$ (63) $T=2{\sqrt {2r\epsilon _{c}-\epsilon _{c}^{2}}}$

Donde ${\textstyle r}$ es el radio de curvatura y ${\textstyle T}$ es el tamaño óptimo a generar en (63). Para el caso de las superficies es conveniente proyectar la derivada segunda en la dirección de la normal. La razón de ello estriba en que nos interesa calcular la curvatura en la dirección perpendicular a la curva y despreciar la curvatura en las direcciones laterales. Con este criterio, una superficie plana con un lado curvo en su mismo plano, tendrá elementos grandes en su interior y éstos se harán sólo pequeños al acercarse a dicha curva. El cálculo quedará de la siguiente manera:

${\vec {D}}_{u}=\left({\frac {\partial _{u}^{2}{\vec {S}}(u,v)}{\partial u^{2}}}\cdot {\vec {N}}\right){\vec {N}}$ (64) $r_{u}={\frac {\left\|{\frac {\partial _{u}{\vec {S}}(u,v)}{\partial u}}\right\|^{3}}{\left\|{\frac {\partial _{u}{\vec {S}}(u,v)}{\partial u}}\times D_{u}\right\|}}$ $T_{u}=2{\sqrt {2r_{u}\epsilon _{c}-\epsilon _{c}^{2}}}$

Siendo ${\textstyle {\vec {N}}}$ la normal a la superficie en (64). El resto de términos son la extensión de las ecuaciones para el caso unidimensional. El cálculo del tamaño óptimo en ${\textstyle v}$ se realiza de manera similar.

Si la superficie se va a mallar de manera no estructurada, se escoge ${\textstyle T=\min(T_{u},T_{v})}$. De esta manera nos aseguramos que se cumple siempre el criterio.

En caso de generación de malla estructurada, se escoge para cada una de las dos direcciones el tamaño calculado. Este forma de mallado permite generar elementos altamente distorsionados en una dirección para superficies con alta curvatura en una sola de sus direcciones principales. Este caso se da típicamente en radios de acuerdo entre superficies. La generación mediante este criterio es especialmente útil para superficies no conformes entre sí, pues, en caso contrario, aparecen los problemas de compatibilización de elementos entre superficies vecinas. Puede comprobarse en el Apartado 4.8.2 cómo aparecen los criterios de incompatibilidad.

En la Figura 35, se puede apreciar el gran ahorro de elementos que se produce en la generación de malla estructurada no conforme de una superficie a otra.

Figura 35: En estas dos mallas se compara la generación por error cordal para malla no estructurada, con la de malla estructurada con superficies independientes.

4.7.2 Corrección de tamaños

Un factor determinante para el éxito en la generación de malla proviene del hecho que los tamaños de elementos no sean excesivamente discordantes entre zonas cercanas. Típicamente, si se contraponen elementos unos contra otros donde la relación de tamaños sea mayor que un factor de 3 ó 4, los algoritmos de puntos cercanos fracasarán en encontrar todas las entidades necesarias para un mallado correcto. Una de las razones proviene del necesario ahorro en la búsqueda, para mejorar la eficiencia computacional. Otro factor es el hecho de que ya no puede crearse una malla bien formada a partir de elementos de tamaño tan discordante.

En el contorno (a) de la Figura 36, en caso de poder generarse una malla, ésta daría elementos de muy mala calidad. El contorno corregido de la Figura (b) proporcionaría una malla correcta.

Figura 36: En la Figura (a) se observa la dificultad de generar una malla en el interior del contorno. En la (b) se corrige tamaños para obtener mejores elementos.

La operación de corrección de tamaños se realiza antes de iniciar la generación de malla y consiste en cambiar internamente los tamaños asignados a las entidades por el usuario.

Un primer paso se produce superficie a superficie, donde cada línea origina un tamaño máximo en función de su longitud. Este tamaño debe traspasarse, adecuadamente incrementado a las líneas vecinas en las superficies a las que pertenecen.

En un segundo paso, para cada entidad que tenga un tamaño menor que el de generación, deben buscarse las entidades suficientemente cercanas y corregir su tamaño asignado. Para ello se escoge un punto definidor de cada entidad, que puede ser su centro, y a partir de él se calculan distancias con el resto de entidades. Dado que es una búsqueda de todos contra todos, es necesario aplicar un método eficiente de búsqueda espacial. El método escogido es el de la búsqueda por Octtree, que será definido en detalle en el Apartado 4.12.1.

El criterio de crecimiento de tamaños con la distancia se basa en una progresión aritmética creciente tomando como base la separación entre centros definidores. Así, si ${\textstyle d}$ es la distancia entre las dos entidades y ${\textstyle S_{0}}$ es el tamaño en la entidad base, el tamaño para la nueva entidad se calculará como:

$S_{1}={\frac {(F-1)d+S_{0}}{F}}$

(65)

donde ${\textstyle F}$ es el factor de crecimiento, que se situará típicamente en valores del orden 1.5 a 2. Debe tenerse en cuenta que este nuevo tamaño sólo se aplicará en caso de ser menor al que la entidad ya tenía previamente asignado.

Por su misma naturaleza, este proceso debe ser iterativo para poder traspasar secuencialmente los tamaños de unas entidades a otras. Habitualmente, con 2 ó 3 iteraciones se consiguen unos buenos resultados.

4.8 GENERACIÓN DE MALLA ESTRUCTURADA

La condición de malla estructurada se impone a las entidades geométricas y requiere que éstas cumplan las siguientes restricciones:

• Toda línea puede tener malla estructurada. En el caso de líneas, no tiene sentido de hablar de malla estructurada, pero se hace la distinción por la manera de definir los tamaños de elemento. En malla no estructurada se usa el esquema general de tamaño de elemento. En cambio, en malla estructurada se define el número de elementos que se desean sobre la línea. Toda línea que sea contorno de una superficie estructurada debe ser definida según este criterio aunque otra superficie de la que sea contorno sea no estructurada.

• Sólo las superficies de cuatro lados o topológicamente rectangulares pueden ser malladas de forma estructurada.

• Los volúmenes deben tener seis caras o superficies de contorno.

Estas restricciones son necesarias para que el algoritmo de generación de mallas estructuradas de cuadriláteros o de hexaedros tenga sentido. Se podrían definir otros algoritmos más sofisticados, especialmente para el caso de triángulos, pero con los existentes se pueden realizar gran variedad de configuraciones.

Figura 37: Mediante el empleo de superficie NURBS de más de 4 lados, se pueden hacer transiciones de tamaño de una superficie a otra.

La restricción más grave viene dada por el hecho de que los lados opuestos de una superficie o volumen deben tener el mismo número de elementos¹. Esto implica que el número de elementos se propaga entre superficies o volúmenes y dificulta la concentración de elementos en zonas determinadas. Puede verse en la Figura 37, que mediante el empleo de NURBS puede llegarse a una definición aceptable de las concentraciones de los elementos.

Cuando se seleccionan algunas superficies o volúmenes, automáticamente se seleccionan todas las líneas que son contorno de estas entidades y el usuario puede asignar número de elementos en cada línea. Cuando ese número queda asignado para una línea, se traspasa a todas las otras líneas que requieran compatibilidad con la primera. Para el caso especial de polilíneas con superficies múltiples, es necesario un algoritmo posterior de comprobación de compatibilidad.

(¹) Esta restricción puede en parte subsanarse con el uso de polilíneas que se describirá en apartados posteriores.

4.8.1 Malla estructurada con polilíneas de superficie múltiple

El caso de polilíneas de superficie múltiple se trata de manera especial y utilizando un algoritmo posterior al de la asignación de datos de malla estructurada. Como en una polilínea de este tipo habrá un número distinto de superficies por un lado que por el otro, no se puede conseguir a priori asegurar que el número de elementos para un lado y para el otro de la línea va a ser el mismo. Para solventar este problema es necesario comparar en cada una de estas polilíneas el número de elementos por un lado y por otro. En caso de ser diferentes, deben compatibilizarse y extender el cambio. El algoritmo desarrollado va actuando iterativamente hasta que compatibiliza todas las superficies estructuradas o comprueba que es imposible (ver la Figura 38). En este último caso dará un mensaje al usuario con indicaciones de que determinada conexión en una polilínea debe arreglarse manualmente.

4.8.2 Malla estructurada en superficies NURBS de más de 4 lados

El mismo algoritmo de unión simple de lados opuestos puede ser aplicable para generar mallas sobre superficies de más de 4 lados, siempre que estas se conserven topológicamente similares a un cuadrilátero. A efectos prácticos esto significa que tengan 4 ángulos pronunciados y el resto de ángulos entre líneas sean nulos o pequeños. Un ejemplo de esta configuración es el de la Figura 37.

Figura 38: Ejemplo de configuración de superficies en el que no puede conseguirse la coherencia en el número de particiones en lados.

Para superficies de este tipo, en el momento de asignar el número de particiones a las líneas no puede calcularse automáticamente el número de particiones para las líneas asociadas. La solución desarrollada consiste en aplicar un algoritmo de comprobación a posteriori de la coherencia en los tamaños. Pueden darse casos, como el de la Figura 38, en el que no puede conseguirse la disposición correcta de particiones.

4.8.3 Generación de los nodos interiores en una superficie estructurada

Todas las superficies que pueden mallarse de manera estructurada deben ser paramétricas y pueden ser trimadas. La función que las define es:

$S\left(u,v\right)\quad {\mbox{con }}\left\{{\begin{array}{l}0\leq u\leq 1\\o\leq v\leq 1\end{array}}\right.$

(66)

Disponemos de los nodos del contorno sobre las líneas en el espacio paramétrico ${\textstyle u,v}$. Si llamamos ${\textstyle {\vec {u}}_{ij}\in R^{2}}$ con ${\textstyle i=0}$ o ${\textstyle j=0}$ a los puntos del contorno, el cálculo de los nodos interiores de la malla estructurada consta de dos fases:

1. El cálculo de ${\textstyle {\vec {u}}_{ij}}$ se realiza a partir de la interpolación Coon explicada en el Apartado 8.9.2. Obtenemos tras ello un vector perteneciente al espacio paramétrico.
2. El nodo en el espacio real se obtiene aplicando la función ${\textstyle {\vec {S}}(u,v)}$.

Este esquema generará elementos con buena relación de aspecto siempre que la superficie sea topológicamente rectangular y no esté excesivamente distorsionada. En la Figura 39, se puede ver un ejemplo de la aplicación de esta técnica.

Figura 39: Generación de malla estructurada en una superficie trimada.

4.8.4 Generación de los nodos interiores en un volumen estructurado

Los volúmenes que pueden mallarse como estructurados deben tener seis caras o superficies de contorno. Cada una de estas superficies de contorno podrá representarse según el algoritmo explicado en el apartado anterior. Por tanto, en el caso de desear la generación de una malla de ${\textstyle M\times {}N\times {}O}$ elementos, el volumen puede representarse como:

$V\left(u=0,v,w\right)=S_{1}(v,w)$

(67)

$V\left(u=1,v,w\right)=S_{2}(v,w)$

(68)

$u_{i}\quad {\hbox{con }}\quad i=1,\dots ,M\quad {\hbox{, }}v_{j}\quad {\hbox{con }}\quad j=1,\dots ,N\quad {\hbox{y }}\quad w_{k}\quad {\hbox{con }}\quad k=1,\dots ,O$

(69)

para todos los contornos¹. El cálculo de un punto interior del volumen se puede realizar mediante varios algoritmos. El más sencillo consiste en la interpolación básica ponderada con la distancia, de los puntos de influencia del contorno y puede expresarse como:

$V_{1}\left(u_{i},v_{j},w_{k}\right)=\left(1-f_{1}\right)\cdot {}S_{1}\left(u_{i},v_{j}\right)+f_{1}\cdot {}S_{2}\left(u_{i},v_{j}\right)\quad {\hbox{con }}f_{1}={\frac {w_{k}}{O}}$ $V_{2}\left(u_{i},v_{j},w_{k}\right)=\left(1-f_{2}\right)\cdot {}S_{3}\left(u_{i},w_{k}\right)+f_{2}\cdot {}S_{4}\left(u_{i},w_{k}\right)\quad {\hbox{con }}f_{2}={\frac {v_{j}}{N}}$ $V_{3}\left(u_{i},v_{j},w_{k}\right)=\left(1-f_{3}\right)\cdot {}S_{5}\left(v_{j},w_{k}\right)+f_{3}\cdot {}S_{6}\left(v_{j},w_{k}\right)\quad {\hbox{con }}f_{3}={\frac {u_{i}}{M}}$ $V\left(u_{i},v_{j},w_{k}\right)={\frac {1}{3}}\sum _{i=1}^{3}V_{i}$

Otra técnica que suele dar mejores resultados consiste en aprovechar la interpolación Coon de manera similar a la empleada en la generación de malla estructurada sobre superficies. En este caso:

${\vec {V}}_{1}(u_{i},v_{j},w_{k})={\vec {S}}(u_{i},v_{j})\quad {\hbox{con }}{\vec {S}}{\hbox{ definida en }}w=w_{k}{\hbox{ }}$

(70)

Siendo ${\textstyle {\vec {S}}}$ el interpolador Coon definido en el Apartado 8.9.2. De manera equivalente se pueden definir ${\textstyle {\vec {V}}_{2}}$ y ${\textstyle {\vec {V}}_{3}}$ sin más que cambiar el plano de interpolación de ${\textstyle w}$ fijado a ${\textstyle u}$ ó ${\textstyle v}$ fijado. Finalmente se puede escoger uno de los ${\textstyle {\vec {V}}_{i}}$ o promediarlos,

$V\left(u_{i},v_{j},w_{k}\right)={\frac {1}{3}}\sum _{i=1}^{3}V_{i}$

El la Figura 40, se describe gráficamente la técnica anterior.

Figura 40: El cálculo de la interpolación volumétrica se realiza mediante la interpolación Coon en algunas de sus secciones transversales.

Para otros métodos más suaves de interpolación de puntos interiores se podrían usar las funciones de forma de un elemento hexaédrico de alto grado (ver [18]) o usar una ecuación diferencial euleriana y resolver el sistema (ver [21]).

(¹) En todo momento deben tenerse en cuenta las orientaciones y disposición de las superficies.

4.8.5 Concentración de elementos en mallas estructuradas

Para algunos tipos de análisis es necesario concentrar elementos en algunas zonas. Típicamente, conviene poner mayor número de elementos en zonas de altos gradientes de la variable de cálculo. Para ello, conviene realizar esta concentración en el espacio paramétrico y actuar únicamente en las líneas. El traspaso hacia las superficies y volúmenes de esta concentración será automático, según los algoritmos de interpolación ya descritos.

Cada línea, por ser estructurada debe tener como definición de tamaño el número de particiones sobre ella. La concentración de elementos la definimos como una progresión geométrica de razón ${\textstyle r}$¹. Con estas condiciones, el posicionamiento de los nodos en el espacio paramétrico será:

$u_{0}$ $=0$ $u_{i}$ $={\frac {f_{n}\cdot u_{i-1}}{u_{n}}}\qquad i=1,\dots ,n$ $f_{i}$ $=r\cdot f_{i-1}$

Figura 41: Diferentes criterios de concentración de elementos para una misma geometría y para el mismo número de particiones.

(¹) La razón $r$ puede tomar un valor positivo o negativo que luego se escala a menor que uno para los valores negativos y mayor que uno para los positivos. El motivo de este escalado es dar facilidad intuitiva al usuario.

4.9 GENERACIÓN DE MALLA EN LÍNEAS

Antes de poder generar malla sobre las superficies y volúmenes, es necesario generar elementos de 2 ó 3 nodos sobre todas las líneas que sean contorno de superficies. Por defecto, estos elementos serán únicamente parte de la malla final si pertenecen a líneas aisladas o si el usuario lo requiere.

Como todos los tipos de línea están definidos de forma paramétrica en el intervalo ${\textstyle \left[0,1\right]}$, se evalua la función ${\textstyle L\left(u\right)}$ con ${\textstyle u\in \left[0,1\right]}$ para diferentes valores del parámetro. Para líneas rectas y arcos, como el parámetro de la curva es el parámetro arco, una vez definidos los intervalos en el espacio del parámetro, quedan bien espaciados los nodos en el espacio real. Para el caso de las líneas NURBS, es necesaria una función de corrección para correlacionar el avance en el espacio del parámetro con el espacio real. Esta función se describirá en apartados posteriores. Las polilíneas son un caso mixto de los anteriores.

4.9.1 Cálculo de los parámetros de espaciamiento de nodos

Para todos los tipos de línea excepto la polilínea, a partir de los tamaños deseados de elemento se calculan las posiciones del parámetro en el espacio ${\textstyle \left[0,1\right]}$. Para ello el algoritmo es: Se dispone de tres tamaños de elementos que están situados en los extremos de la línea y en su centro. Estos tamaños provendrán del tamaño asignado por el usuario o por algún algoritmo de asignación o por defecto. En el caso de mallas estructuradas, como el dato es número de elementos por línea, se calculan los tamaños unitarios:

$T_{1}=T_{2}=T_{3}={\frac {1.0}{NE}}$

(71)

siendo ${\textstyle T_{i}}$ los tamaños en el espacio ${\textstyle \left[0,1\right]}$ y ${\textstyle NE}$ el número de elementos sobre la línea. En malla no estructurada, cada tamaño en el espacio debe convertirse a ese mismo espacio:

$T_{i}={\frac {S_{i}}{L}}$

(72)

siendo ${\textstyle S_{i}}$ el tamaño deseado y ${\textstyle L}$ la longitud de la línea. Para posicionar los parámetros se parte de una configuración de espaciamiento constante y separación del mínimo de los tres valores ${\textstyle T_{0}=\min \left(T_{1},T_{2},T_{3}\right)}$. A partir de ella se va iterando partiendo de la base de que se tiene una función ${\textstyle S\left(t\right)}$ tal que nos da el espaciamiento óptimo para cada valor de ${\textstyle t}$

$t_{i}=t_{i-1}+S\left({\frac {t_{i}+t_{i-1}}{2}}\right)$

(73)

Figura 42: Cambio del tamaño óptimo entre $t=0$ y $t=0.5$ para diferentes factores de transición.

Finalmente se corrigen los parámetros con un factor

$F={\frac {1}{t_{NE}+S\left({\frac {1+t_{NE}}{2}}\right)}}$

(74)

$t_{i}=t_{i}\cdot {}F$

(75)

La definición de la función ${\textstyle S\left(t\right)}$ viene dada por el hecho de que tiene que haber una transición adecuada entre tamaños pequeños y grandes y que esta transición viene regulada por el factor de transición definido por el usuario. Su formulación viene dada por dos parábolas, una entre ${\textstyle t=0}$ y ${\textstyle t=0.5}$ y la otra entre ${\textstyle t=0.5}$ y ${\textstyle t=1.0}$ Si ${\textstyle T_{1}}$ es el tamaño deseado en ${\textstyle t=0}$, ${\textstyle T_{2}}$ es el tamaño en ${\textstyle t=0.5}$ y ${\textstyle SF}$ es el factor de transición definido en ${\textstyle \left[0,1\right]}$ entonces:

$T_{tg}=T_{\min }+SF\cdot {}\left(T_{\max }-T_{\min }\right)\quad {\hbox{con }}\left\{{\begin{array}{l}T_{\min }=\min \left(T_{1},T_{2}\right)\\T_{\max }=\max \left(T_{1},T_{2}\right)\end{array}}\right.$

(76)

$S\left(t\right)=T_{1}\left(1-2t\right)^{2}+2(2t)(1-2t)T_{tg}+T_{2}\left(2t\right)^{2}\quad {\mbox{con }}t\in \left[0,0.5\right]$

(77)

Un ejemplo de esta parábola puede verse en la Figura 42.

La parábola en ${\textstyle t\in \left[0.5,1\right]}$ es la misma que la anterior, adecuadamente trasladada.

4.9.2 Corrección para líneas NURBS

Dado que todas las operaciones anteriores se han realizado en el espacio del parámetro ${\textstyle t\in \left[0,1\right]}$ y dado que lo que en realidad se desea es mantener las mismas relaciones en el espacio real, sólo se obtendrá un resultado correcto si el parámetro de la curva es el parámetro arco. Esto se cumple para las líneas rectas y para los arcos pero no para las NURBS (en la Figura 43, se puede apreciar un ejemplo típico de mala parametrización). En este último caso será necesario hacer una corrección para que las relaciones en un espacio se mantengan en el otro. El algoritmo es el siguiente: Dada una lista de parámetros ${\textstyle t_{1},\dots ,t_{n}\in \left[0,1\right]}$ se desea obtener otra lista ${\textstyle t_{1}^{\prime },\dots ,t_{n}^{\prime }\in \left[0,1\right]}$ tal que las relaciones de distancia ${\textstyle {\frac {|t_{i}-t_{j}|}{1}}}$ se mantengan en el espacio ${\textstyle {\frac {|{\vec {L}}\left(t_{i}\right)-{\vec {L}}\left(t_{j}\right)|}{L_{tot}}}}$. El proceso consiste en calcular cada ${\textstyle t_{i}^{\prime }}$ como:

$t_{i}^{\prime }=t_{i-1}^{\prime }+I\left(t_{i-1}^{\prime },\Delta t_{i-1,i}\right)$

(78)

y corregir la nueva lista obtenida para que cubra el espacio ${\textstyle \left[0,1\right]}$:

$t_{i}=t_{i}\cdot {}{\frac {1}{t_{NE}+S\left({\frac {1+t_{NE}}{2}}\right)}}$

(79)

La función de incremento óptimo ${\textstyle I\left(t_{i-1}^{\prime },\Delta t_{i-1,i}\right)}$ viene dada por el concepto de que:

$t_{i}^{\prime }=t_{i-1}^{\prime }+{\frac {\int _{t_{i-1}^{\prime }}^{t_{i}^{\prime }}\left|{\frac {dL}{dt}}\right|dt}{L_{tot}}}$

(80)

Figura 43: Para líneas rectas y para arcos, la parametrización de la línea viene dada por el parámetro arco. En NURBS, puede no ser así.

Como el límite superior de la integral no se conoce, se va iterando a partir del incremento de partida y se calcula la integral con una cuadratura de Gauss-Legendre de número de puntos dependiente de la iteración. Conocido ${\textstyle t_{i-1}^{\prime }}$ y ${\textstyle \Delta t_{0},}$ se calcula ${\textstyle t_{i}^{\prime }=t_{i-1}^{\prime }+\Delta t_{0}}$ y:

$\Delta L=\int _{t_{i-1}^{\prime }}^{t_{i}^{\prime }}\left|{\frac {dL}{dt}}\right|dt=\sum _{i=1}^{n_{p}}\omega _{i}\left.{\frac {dL}{dt}}\right|_{t_{GL}}$

(81)

${\mbox{si }}\left|\Delta L-\Delta t_{i}L_{tot}\right|<\varepsilon \quad {\mbox{se para. Si no }}$

(82)

$\Delta t_{i+1}={\frac {\Delta L}{L_{tot}}}$

(83)

En todo este proceso debe haber controles para casos extremos de ${\textstyle \Delta t}$ muy altos y curvas muy mal parametrizadas.

La corrección dada en (79), puede dar lugar a grandes errores para líneas mal parametrizadas. La razón es que estamos multiplicando cada ${\textstyle t_{i}}$ con un valor constante para todos ellos. Por tanto desplazamos por igual todos los valores, tanto por zonas de módulo de la derivada pequeña como por zonas de derivada grande. Una mejor corrección es:

$CF_{0}={\frac {1}{t_{NE}+S\left({\frac {1+t_{NE}}{2}}\right)}}$

(84)

$CF=1.0+\left(CF_{0}-1.0\right)\cdot {}{\frac {t_{i}^{\prime }-t_{i-1}^{\prime }}{t_{i}-t_{i-1}}}$

(85)

$t_{i}^{final}=t_{i-1}+CF\cdot {}\left(t_{i}^{\prime }-t_{i-1}^{\prime }\right)$

(86)

Este valor debe tener los oportunos controles de corrección de valores por encima del parámetro 1.0

4.9.3 Generación para polilíneas

Las polilíneas deben tratarse de manera especial ya que son contenedoras de otros tipos de líneas y tienen toda una casuística sobre si los puntos interiores de la malla deben convertirse en nodos o no. Las posibilidades que se han implementado son:

• Ningún punto forma parte de la malla. En este caso la línea puede considerarse como cualquier otra.

• Algunos o todos los puntos forman parte de la malla. Aquí, para definir los tamaños de mallado, se escoge el conjunto de líneas entre dos puntos que deben ser mallados y como tamaños extremos, se escoge los de estos dos puntos. El tamaño intermedio no puede ser el de la línea porque no permitiría concentrar en diferentes segmentos. Por tanto, el tamaño en el centro de la línea se calcula a partir de los dos tamaños extremos. A partir de aquí, el proceso es el mismo que para una línea normal con el simple escalado y translación de los parámetros.

4.10 GENERACIÓN DE MALLA EN SUPERFICIES. INTRODUCCIÓN

En general, al hablar de generación de malla en superficies, nos referimos a superficies en el espacio. Las superficies en el plano ${\textstyle z=0}$, o usualmente consideradas como 2D, serán tratadas como un caso particular de las anteriores.

Existen varias posibilidades para realizar esta generación:

• Transformar la superficie a una superficie equivalente en un plano y realizar una generación 2D en él.
• Generar la malla directamente en el espacio.
• Generar la malla en el espacio paramétrico con la métrica del espacio real.

Los dos primeros métodos han sido desarrollados en esta tesis y cada cual contiene una serie de ventajas e inconvenientes. En la descripción que sigue, se va a empezar con el método de generación en el plano. En el se describirá detalladamente el método de generación así como las técnicas para realizar la transformación directa y la inversa. En la descripción de la generación en el espacio se comentará únicamente las diferencias con el método anterior. Finalmente se discutirán las ventajas e inconvenientes de cada método.

4.11 GENERACIÓN DE MALLA EN EL PLANO

Seguidamente se describirán las técnicas de generación en el caso de transformación previa al plano. De esta manera se consigue reducir un problema tridimensional a uno bidimensional¹.

(¹) Algunas veces también llamado $2{\frac {1}{2}}D$ por ser superficies en el espacio 3D.

4.11.1 Transformación del espacio paramétrico al espacio equivalente

El paso previo en este método a la generación de la malla propiamente dicha es la transformación de los elementos de contorno o malla sobre línea, a un espacio equivalente que esté contenido en un plano y que tenga la forma adecuada para minimizar la distorsión.

En el caso de superficies planas, la transformación consiste simplemente en una rotación tal que la normal al plano se convierta en ortogonal al plano ${\textstyle z=0}$. Se desprecia la coordenada z, se genera la malla y se aplica la operación inversa.

En superficies Coon o NURBS, al ser paramétricas en el espacio ${\textstyle \left[0,1\right]\times {}\left[0,1\right]}$, la opción más sencilla sería generar malla en el espacio paramétrico. La generación directa en este espacio y posterior transformación al espacio real daría lugar a elementos con fuertes distorsiones. Para subsanar este problema hay básicamente dos alternativas:

• Generar malla en el plano con la métrica del espacio real. Básicamente se trata de que todas las operaciones que requieran una distancia o un ángulo, en lugar de realizarlas directamente se efectuan en la superficie real mediante unos operadores de transformación. Ver [8] para una descripción más detallada de este esquema.

• Convertir el espacio paramétrico, que es un cuadrado, en una figura geométrica más similar a la superficie real de manera que haya un parecido mayor de ángulos y de relaciones de distancias.

• Un modelo mixto entre ambos.

El esquema escogido ha sido el de convertir el espacio paramétrico a una forma geométrica más adecuada. La ventaja radica en que se evitan en el generador las dificultades adicionales de tratar con estos operadores de transformación. Se separa así la problemática de la generación con la de la transformación. Con esta técnica se consigue muy poca distorsión de los elementos para formas muy variadas de las superficies. Otra alternativa de posible consideración en el futuro sería una solución mixta entre ambas posibilidades.

La definición de esta nueva forma geométrica viene dada por asociación a un elemento cuadrilátero de los usados habitualmente en elementos finitos. Puede usarse tanto uno de cuatro nodos como uno de ocho. En algunos casos degenerados, este cuadrilátero se convertirá en triángulo para adaptarse mejor a la forma geométrica real. Este cuadrilátero se define de tal manera que se minimice la distorsión de los elementos resultantes. Para ello el algoritmo es el siguiente:

Figura 44: Relación entre la superficie real y el espacio donde se realiza la generación de malla bidimensional.

Dados los 4 lados de la superficie Coon o NURBS ${\textstyle l_{i}}$, y los 4 ángulos ${\textstyle \alpha _{i}}$, se escalan las longitudes obteniéndose ${\textstyle {\bar {l}}_{i}={\frac {l_{i}}{l_{o}}}}$. Los ángulos se escalan para sumar ${\textstyle 2\pi }$ radianes. La matriz de definición del cuadrilátero será:

$G=\left[{\begin{array}{llll}0.0&1.0&1.0-{\bar {l}}_{2}\cos \left(\alpha _{1}\right)&{\bar {l}}_{4}\cos \left(\alpha _{4}\right)\\0.0&0.0&{\bar {l}}_{2}\sin \left(\alpha _{1}\right)&{\bar {l}}_{4}\sin \left(\alpha _{4}\right)\end{array}}\right]$

(87)

Para obtener los nodos de centro de cara del cuadrilátero de 8 nodos, se calculan primero como el punto medio de cada uno de sus dos extremos. Después se les aplica un factor de corrección en la dirección normal a la cara. Ese factor se calcula como:

${\vec {V}}={\vec {S}}\left(u=0.5\right)-{\vec {S}}\left(u=0.0\right)-{\frac {{\vec {S}}\left(u=1.0\right)-{\vec {S}}\left(u=0.0\right)}{2.0}}$

(88)

$f={\frac {{\vec {V}}\cdot {}{\frac {\partial S}{\partial v}}\left|_{u=0.5}\right.}{\left|{\vec {V}}\right|\left|{\frac {\partial S}{\partial v}}\left|_{u=0.5}\right.\right|}}$

(89)

En la Figura 44 se puede ver la relación de formas entre la superficie real y el espacio calculado.

Para generar la malla, se obtiene la posición de los nodos de las líneas de contorno en el espacio paramétrico de la superficie y se les aplica la transformación definida por ${\textstyle G}$. Se genera la malla en el plano y se realiza la transformación inversa. Este paso del espacio del cuadrilátero al paramétrico se realiza con un método iterativo como un Newton-Rapson ya que se conoce la matriz jacobiana de la transformación. El paso del espacio paramétrico al real es directo mediante la aplicación de la función que define la superficie.

4.12 ALGUNAS ESTRUCTURAS DE DATOS USADAS EN LOS ALGORITMOS DE GENERACIÓN

Previamente a la descripción de la generación de malla 2D y 3D, se explicarán algunas estructuras de datos usadas en este proceso. Al hablar de estructuras de datos, se entiende que llevan asociados un conjunto de algoritmos para su uso.

4.12.1 Algoritmos de búsqueda de puntos cercanos: Quadtree y Octtree

Dado un conjunto de puntos en el plano o en el espacio, el objetivo es encontrar todos los cercanos a un punto dado. Una de las hipótesis de partida es evitar la comparación directa de distancias entre todos ellos con el punto base. O sea, dado que el algoritmo básico es de orden ${\textstyle O^{2}}$, interesa encontrar otro algoritmo que sea de orden ${\textstyle n\log(n)}$.

Para minimizar el número de comparaciones, se recurre a una estructura de tipo Quadtree en el plano y de tipo Octtree¹ en el espacio. Se va a describir el algoritmo para el caso en el plano. El caso 3D es una extensión directa del caso bidimensional.

La estructura de datos Quadtree se basa en subdividir el plano en un conjunto de regiones para acotar la búsqueda sólo en las regiones de interés. La definición de las subregiones se va realizando a medida que se introducen los puntos. El proceso es el siguiente:

Figura 45: Estructura recursiva del Quadtree. En la figura de la izquierda se observa la forma geométrica. En la de la derecha, la forma de árbol resultante.

1. Se parte de una única región que es un rectángulo que cubre todo el dominio. Para ello debe conocerse con antelación su tamaño. Normalmente se itera previamente sobre todos los puntos para calcular la región mínima que los contiene a todos. Existen técnicas para aumentar dinámicamente la región de partida que serán descritas más adelante.
2. Se van introduciendo puntos en la región definida. Estos puntos se almacenan en esta región hasta llegar a un máximo permitido por región (este máximo es arbitrario y pueden usarse 4 puntos por región en los Quadtree y 8 en los Octtree).
3. Cuando la región está llena, ésta se subdivide en 4 subregiones que se crean dividiendo el rectángulo original en cuatro rectángulos iguales. Tanto los puntos que ya estaban en la región, como el nuevo punto entrado, se introducen en la subregión a la que pertenecen en función de sus coordenadas. Los sucesivos puntos que se vayan introduciendo, se alojarán en su subregión correspondiente, en el nivel más bajo posible. Por tanto, se entiende que una región determinada aloja o puntos o subregiones pero no ambas cosas a la vez.
4. Se procede de manera recursiva a medida que se van entrando todos los puntos y se van creando todas las subregiones que sean necesarias. De esta manera se llega a una estructura de árbol.
5. Cuando todos los puntos están introducidos, el Quadtree está creado y lleno. A partir de aquí, se puede proceder a la búsqueda.

La forma típica de un Quadtree puede observarse en la Figura 45. El funcionamiento del Octtree es equivalente al descrito, cambiando los rectángulos por hexaedros en el espacio. Un ejemplo de una estructura de este tipo puede observarse en la Figura 46.

Figura 46: Esquema de la forma geométrica típica de un Octtree.

Mediante una estructura de este tipo se pueden realizar varias operaciones. La más común de ellas es, dado un punto de partida, encontrar un conjunto de puntos cercanos. El criterio de cercanía puede ser tanto una distancia máxima como un número óptimo o máximo de puntos. También puede ser un criterio que contemple las dos opciones anteriores.

4.12.1.1 Puntos cercanos al punto dado según un radio máximo

En algunos casos puede interesar encontrar todos los puntos cuya distancia a un punto de referencia sea menor a un radio dado. Para ello se define una función recursiva que empieza en la raíz y que actuará de manera diferente según la rama sea contenedora de hijos o contenga directamente puntos. Si llamamos ${\textstyle H_{i1}}$ al hijo con ${\textstyle x_{i}}$ menor, ${\textstyle H_{i2}}$ al hijo con ${\textstyle x_{i}}$ mayor y ${\textstyle c_{i}}$ al centro del contenedor según ${\textstyle x_{i}}$, el proceso será:

$f:\left\{{\begin{array}{rl}{\mbox{contiene hijos: }}&\left|{\begin{array}{ll}{\mbox{si }}P_{i}-R\leq c_{i}&{\mbox{se entra en }}H_{i1}{\mbox{ }}\\{\mbox{si }}P_{i}+R\geq c_{i}&{\mbox{se entra en }}H_{i2}{\mbox{ }}\end{array}}\right.\\{\mbox{contiene puntos: }}&{\mbox{se acepta }}{\vec {P}}_{I}{\mbox{ si }}\|{\vec {P}}_{I}-{\vec {P}}\|\leq R\end{array}}\right.$

Siendo ${\textstyle P}$ el punto de referencia, ${\textstyle {\vec {P}}_{I}}$ el punto interior a comprobar y ${\textstyle R}$ el radio máximo aceptado.

4.12.1.2 Numero fijo de puntos más cercanos a uno dado

Aunque sería posible definir una función que diera un número exacto de puntos alrededor de uno dado y que cumplieran el criterio de ser los más cercanos a él, esta función sería muy ineficiente en tiempo de computación. Por otra parte, en muchos casos deseamos un número aproximado de puntos cercanos a uno dado y que estén contenidos en un radio máximo.

4.12.1.3 Número aproximado de puntos cercanos a uno dado

Para la generación de malla por el método del avance frontal interesa encontrar los puntos cercanos a la cara activa del avance. Para ello tenemos una primera restricción que es que sólo es necesario considerar el semiespacio definido por la cara y su orientación. O sea que si ${\textstyle {\vec {P}}}$ es el centro de la cara y ${\textstyle {\vec {N}}}$ la normal a la cara, los puntos ${\textstyle {\vec {P}}_{I}}$ a encontrar deben cumplir que ${\textstyle ({\vec {P}}_{I}-{\vec {P}})\cdot {\vec {N}}\geq \varepsilon }$.

Por otro lado, la noción de puntos cercanos se puede considerar relativa. Debido a los cambios bruscos de tamaño entre elementos de la malla, es difícil calcular un radio de cercanía en el que se encuentren todos los puntos de interés. Hay que tener en cuenta que encontrar un número demasiado pequeño de puntos supondrá el fracaso del método. En cambio, una elección conservadora en exceso, provocará un incremento del coste de ejecución que se convertirá en inabordable fácilmente. En la Figura 47, se comprueba la diferencia de tamaño del área de búsqueda para dos configuraciones del frente similares.

Figura 47: En estas figuras se puede observar la diferencia de tamaño del área de búsqueda para una misma cara del frente.

La solución más efectiva consiste en definir un número óptimo aproximado de puntos que hay que encontrar alrededor de uno dado. Este número se calcula en función de los puntos que hay activos en el frente en una situación extrema y que puedan afectar a la creación del nuevo elemento.

El proceso de búsqueda de este número óptimo de puntos es el siguiente: dado un radio máximo ${\textstyle R_{\mathrm {max} }}$, un número óptimo de puntos ${\textstyle N}$, la normal al elemento ${\textstyle {\vec {N}}}$ y el coeficiente ${\textstyle D}$ del plano del elemento, definido tal que un punto ${\textstyle {\vec {P}}}$ estará en el semiespacio de búsqueda si:

${\vec {P}}\cdot {\vec {N}}+D\geq \varepsilon$ (90) entonces el algoritmo de búsqueda será tal que dado un octante hijo $H_{i}$, se descartará si se cumplen las dos condiciones siguientes: ${\vec {P}}_{j}+R_{\mathrm {max} }<H_{i,j,min}$ (91) ${\vec {P}}_{j}-R_{\mathrm {max} }>H_{i,j,max}$ $j={x,y,z}$ Dado un octante padre $H_{p}$, se pueden dar tres casos: $\left\{{\begin{array}{l}H_{p}\quad {\hbox{contiene puntos. (a) }}\\H_{p}\quad {\hbox{contiene un solo octante válido. (b) }}\\H_{p}\quad {\hbox{contiene varios octantes válidos. (c) }}\end{array}}\right.$

(94)

En el caso (92.a), se entrará en el algoritmo 95 con el radio máximo y sin radio mínimo. En el caso (92.b), se procederá recursivamente con la misma función 92 en ese octante. El caso (92.c) requiere un procedimiento iterativo sobre el algoritmo 95.

Esta iteración consta de la entrada en el nuevo algoritmo con un radio máximo menor al dado. Si se considera un máximo de ${\textstyle r}$ iteraciones, se entrará inicialmente con un radio máximo de ${\textstyle R_{max}^{0}=R_{max}/r}$. Si el número de puntos así obtenidos es igual o superior al óptimo, se finaliza la función. En caso contrario, se vuelve a entrar con ${\textstyle R_{max}^{1}=2R_{max}/r}$ y ${\textstyle R_{min}^{1}=R_{max}^{0}}$. El uso del radio mínimo viene dado para evitar la repetición de puntos en las sucesivas iteraciones.

El siguiente criterio a aplicar sobre los octantes hijos ${\textstyle H_{i}}$, dados los datos del algoritmo anterior y un dato adicional ${\textstyle R_{min}}$ es:

$\left\{{\begin{array}{l}H_{i}\quad {\hbox{contiene puntos. (a) }}\\H_{i}\quad {\hbox{es válido según 91 y }}R_{min}{\hbox{ (b) }}\\H_{i}\quad {\hbox{es válido según 91 y }}R_{max}{\hbox{ (c) }}\\H_{i}\quad {\hbox{no es válido según 91. (d) }}\end{array}}\right.$

(98)

En el caso (95.a), se aceptarían todos los puntos que cumpliesen 90. En el caso (95.b), no se entraría en el octante, pues ya se ha usado en la iteración anterior. En el (95.c), se ejecutaría recursivamente la misma función. El último caso, el (95.d), implicaría descartar ese octante.

Al finalizar estas operaciones, se habrá obtenido un número igual o mayor de puntos al óptimo, siempre que el radio máximo lo permita.

La razón de requerir un algoritmo tan complejo para optimizar el número de puntos obtenido, es que una apreciación por defecto provoca un error² irrecuperable en la generación de la malla. Todos los algoritmos de la generación como la ordenación de puntos, cálculo de intersecciones y otros van referidos al número de puntos cercanos, entonces, una apreciación en exceso de ese número, conllevará aumentos en el tiempo de computación directamente proporcionales a ese incremento de número de puntos.

4.12.1.4 Incremento del tamaño del Octtree

En algunos casos se conoce a priori el tamaño máximo de la región donde se encuentran todos los puntos. Esto sería el caso típico de generación de malla bi o tridimensional donde se parte del contorno externo. Este contorno define perfectamente una región convexa que contendrá a todos los puntos. En estos casos no es necesario incrementar el tamaño del Octtree. Incluso hay que evitar hacerlo para poder detectar rápidamente errores en la generación de malla.

En otros casos, como por ejemplo en la generación de malla de superficie en el espacio, el contorno no define esa región convexa contenedora. Entonces es necesario usar un mecanismo de incremento de tamaño del Octtree. El mecanismo es el siguiente: dado un punto ${\textstyle {\vec {P}}}$ que no pertenece a la región mayor del Octtree, se comprueba a que lado de la estructura reside. Entonces se crea un contenedor tal que uno de sus octantes sea la región actual. El resto de las regiones se escogen de manera que la estructura crezca hacia el lado en que está el punto. Si la nueva estructura aún no contiene al punto, se procede recursivamente hasta obtener una que si lo contenga. En la Figura 48, se puede ver un caso típico de recrecimiento en el que sólo es necesario un nivel de aumento.

Figura 48: Recrecimiento de un Quadtree al insertar un punto que no pertenece al interior de la estructura existente.

4.12.1.5 Control de puntos coincidentes

La estructura anteriormente definida funcionará perfectamente si todos los puntos introducidos son distintos entre si. En el caso de que haya varios puntos coincidentes o casi coincidentes³, pueden aparecer problemas. Concretamente, si el número de puntos coincidentes es superior al número de puntos que puede contener una región, los algoritmos recursivos irán iterando hasta el infinito.

Para solventar este problema se define una tolerancia, que puede depender del tamaño del Octtree o puede estar fijada por otros criterios. Cuando se va descendiendo en la estructura y se llega a un octante de dimensión menor que la tolerancia, para todos los puntos que desciendan por ese octante se deja de tener en cuenta el criterio de posicionamiento según coordenadas y se van colocando los puntos en las subestructuras con el único criterio de llenarlas equitativamente. En los algoritmos de búsqueda de puntos cercanos, deberá controlarse esta eventualidad y empezar una búsqueda total en los octantes que entren dentro de esta categoría.

El inconveniente de este acercamiento es que la búsqueda se convierte en orden ${\textstyle O^{2}}$ para los puntos pertenecientes a estos octantes especiales. Raramente pasa a ser un problema real debido a que los puntos que cumplen este criterio suelen ser una cantidad muy baja respecto del total.

(¹) La traducción de Quadtree y Octtree podría ser árbol cuadrante y árbol octante.

(²) Es posible en algunos casos raros recuperarse de ese error.

(³) Es habitual tener puntos coincidentes en la importación o en la transformación de una geometría. También en algunas situaciones especiales de generación de malla.

4.12.2 Estructura de árbol binario

Dado un conjunto de datos que pueden ordenarse mediante una magnitud escalar, el objetivo es mantener una lista eficiente de los elementos de menor a mayor, en la que puedan insertarse y extraerse elementos con facilidad.

Este tipo de lista va a usarse para almacenar las caras del frente de avance en la generación de malla (ver el Apartado 4.14). El proceso consiste en que se va introduciendo nuevas caras de tamaño arbitrario y se borran otras ya usadas. Al mismo tiempo interesa poder acceder siempre a la cara menor, que es la que usaremos para proseguir en el avance.

La característica principal a tener en cuenta es que, en general, habrá muchas caras que tendrán tamaños muy parecidos o iguales entre ellas.

Para la creación de esta lista se readapta el algoritmo de Octtree descrito en el Apartado 4.12.1, a una sola dimensión. La diferencia principal estriba en que si el número de caras iguales es muy grande, el algoritmo Octtree se convierte en orden ${\textstyle O^{2}}$. Es necesario hacer una readaptación que consta de que cada elemento de la lista almacena un puntero a su rama para poder extraerlo de la estructura rápidamente. Como las caras se extraen en el frente siempre de más pequeña a más grande, únicamente optimizaremos el algoritmo de extracción para la cara más pequeña. Por otra parte, la inserción de caras siempre será eficiente porque se seguirá el procedimiento descrito en el Apartado 4.12.1.5.

Esta estructura de árbol binario no debe confundirse con otras parecidas descritas en [22], que consisten en una alternativa al Octtree tridimensional.

4.12.3 Estructuras de nodos, elementos y caras

Para conseguir eficiencia de computación, las búsquedas de entidades debe ser eficiente. El problema de la búsqueda de puntos queda resuelto mediante la utilización de los Octtree ya descritos. Es conveniente, pues, encontrar técnicas de localización de elementos y de caras.

Una primera posibilidad sería la creación de estructuras similares a las ya descritas para la búsqueda de puntos. En este caso, habría que triplicar el coste de búsquedas por tres.

Otra opción más conveniente, que es la que se ha usado, consiste en sacrificar una cierta cantidad de memoria a cambio de disminuir el tiempo de computación. En concreto, mantener apuntadores en algunas de las estructuras que den acceso inmediato a las estructuras próximas. Las conexiones entre entidades serán las de la Tabla 3.

Tabla. 3 Relaciones de conexión entre las entidades.

	=
= Nodos	=
= Elementos	=
= Caras
Nodos	-	Cada nodo apunta a uno de sus elementos.	Cada nodo mantiene todas sus caras activas del frente.
Elementos	Los elementos mantienen una lista de sus conectividades.	Mantiene puntero a todos los elementos que comparten una cara con él.	Apunta a las caras activas del frente que parten de él.
Caras	Las caras mantienen una lista de sus conectividades.	Se apunta al elemento ya generado.	-

Algunos de los algoritmos donde se usará esta información son:

• Al crear un nuevo elemento en el frente, una vez encontrados los puntos cercanos, necesitamos las caras cercanas para el cálculo de intersecciones. Mediante los punteros de nodos a caras activas se obtienen con costo computacional prácticamente nulo.
• En el nuevo elemento creado, se debe incluir el puntero a elementos vecinos. Éste se obtiene a partir de los punteros de caras hacia elementos.
• Al usar la técnica de intercambio de diagonales en el alisado (ver el Apartado 4.16.3.1), se obtiene el elemento contiguo mediante el puntero a los elementos vecinos.
• Uno del los algoritmos del alisado, consiste en el movimiento de nodos. Para poderlo realizar se necesitan los elementos y nodos que tiene alrededor. Éstos se consiguen de la siguiente manera:

1. A partir de un nodo obtenemos uno de los elementos que lo contienen.
2. Los otros elementos provienen de una búsqueda aprovechando la información de elementos contiguos¹.
3. Los nodos se obtienen directamente de las conectividades de los elementos.

(¹) Un nodo podría mantener punteros a todos los elementos que le contienen. No se ha escogido esta opción por la gran cantidad de información a almacenar.

4.12.4 Cubos MinMax

Otra estructura destinada a incrementar la eficiencia de proceso son los cubos MinMax. Éstos se definen como el cubo envolvente de una entidad geométrica. Dada una arista o una cara de un frente en generación de malla tetraédrica, el cubo MinMax se calcula a partir de las coordenadas mínimas y máximas de cada uno de los nodos que las forman.

Los cubos MinMax se usan como primer discriminador en el cálculo de intersecciones entre caras y aristas. En caso de que los respectivos cubos no intersecten, que es una comprobación muy eficiente, no es necesario realizar el cálculo convencional de intersección cara arista. En la Figura 49, se pueden ver ejemplos de ambas posibilidades.

Figura 49: En la Figura (A), se puede asegurar que la cara y la arista no intersectan. En la (B), no se puede decidir la intersección a partir del cubo MinMax y deben usarse los algoritmos de intersección cara-arista.

4.13 MÉTODOS DE GENERACIÓN DE MALLA TRIANGULAR EN 2D Y TETRAÉDRICA EN 3D

Dado que los los principios básicos de generación de malla en 2D para triángulos y en 3D para tetraedros son muy similares, se hará una descripción única para ambos, matizándose las diferencias cuando sea necesario.

Los métodos de generación de malla más populares son:

• Generación frontal (advancing front).

• Delaunay (o Boronoi).

• Generación por particiones, por octtree u otros.

En esta tesis se han desarrollado los dos primeros y se explican con detalle en apartados posteriores. Para más información sobre las otras metodologías de generación de malla ver [21].

4.14 GENERACIÓN FRONTAL

El método de la generación frontal apareció para solventar la necesidad de calcular geometrías lo suficientemente complejas como para no poder discretizarse con malla estructurada. Los primeros estudios para el desarrollo y mejora de esta metodología han sido llevados a cabo por Lohner, Peraire y Peiró (ver [3] y [4]).

La generación frontal (advancing front) se basa en un algoritmo que va creando los elementos uno a uno hasta que cubren todo el interior del objeto a mallar. Este proceso de creación puede esquematizarse en los siguientes pasos:

1. El punto de partida es una malla de contorno. En el caso de generación de malla en el plano, esta malla de contorno está formada por elementos lineales de dos nodos. Para mallas tridimensionales, el contorno está formado por una malla triangular. En el inicio, a esta malla se la llama frente activo. Un ejemplo de este frente puede verse en la Figura 50.
2. El avance debe empezar en un de las caras del frente activo. Se escoge la más pequeña de ellas porque optimiza mejor el futuro tamaño de los elementos.
3. A partir de esa cara se calcula la posición óptima de un punto tal que unido a ella forme un triángulo o tetraedro de tamaño escogido. Este tamaño es una información de entrada para el método y puede venir dado por una malla de soporte (background mesh¹) o puede venir dado en relación a las entidades geométricas que rodean al dominio (ver Apartados 4.6 y 4.7).
4. A la lista de posibles puntos de formación del nuevo elemento, deben añadirse los puntos del frente activo más cercanos al calculado. Esta lista debe ordenarse según un criterio de distancias al punto creado². Todo punto del frente cuya distancia, ${\textstyle D}$ al punto calculado sea menor que ${\textstyle D<f(h)}$ siendo ${\textstyle f(h)}$ una función que depende de la altura del elemento calculado, se considera más óptimo que este punto.
5. Para cada uno de los puntos de la lista ordenada, se calcula si el elemento formado con este punto cumple con una serie de criterios mínimos de calidad, como ángulo mínimo aceptable, etc. También se comprueba si sus caras intersectan con cualquiera de las caras del frente. El primer punto de la lista que cumpla con todos estos criterios queda aceptado y se crea un nuevo elemento.
6. Se actualiza el frente quitando la cara, o caras usadas y añadiendo las nuevas caras que se han creado.
7. Mientras queden caras en el frente, se repite el proceso desde el paso 2${\textstyle ^{o}}$.

Figura 50: Frente activo en una generación de malla de tetraedros mediante el método del avance frontal.

A lo largo de este proceso se va desarrollando una problemática que se va a describir en lo sucesivo. También se describirán las distintas técnicas usadas para resolver estos problemas.

(¹) Malla que cubre el dominio a mallar y que se usa para definir una magnitud, como el tamaño deseado u otras, en cada uno de sus nodos.

(²) En realidad, el punto de referencia para calcular las distancias es ${\vec {P}}_{c}=(1-\alpha ){\vec {C}}_{c}+\alpha {\vec {P}}$, siendo ${\vec {C}}_{c}$ el centro de la cara y $\alpha$ un parámetro ponderador del orden de $\alpha \cong {0.8}$.

4.14.1 Incremento del tamaño de los elementos

A partir de la cara del frente desde la cual se realiza el avance, un dato que hay que precisar es el tamaño óptimo del nuevo elemento¹ a partir de información de tamaños deseados asignada a las entidades geométricas. A efectos prácticos, eso significa que se dispone de información de tamaños en el contorno del dominio a mallar y también se dispone de un tamaño genérico en el interior del dominio. Otras técnicas de cálculo de tamaños óptimos pueden encontrarse en [23].

La idea básica consiste en que el tamaño óptimo se transmite de unos elementos a otros. Resulta interesante fijar una desviación máxima de este parámetro en relación al elemento del frente desde el cual se parte. Esta limitación se define de la siguiente manera:

${\frac {h}{\alpha }}<h_{\mathrm {opt} }<h\cdot \alpha$

(99)

$\alpha =0.66\cdot {}F_{s}^{2}-0.19\cdot {}F_{s}+1.03$

(100)

Siendo ${\textstyle h}$ el tamaño del lado escogido del frente, ${\textstyle h_{\mathrm {opt} }}$ el nuevo tamaño deseado, que puede provenir de una malla soporte o de una definición de tamaños en entidades geométricas, y ${\textstyle F_{s}}$ una factor adimensional con ${\textstyle F_{s}\in \left[0,1\right]}$, que viene dado por el usuario y que representa un incremento lento o rápido del tamaño de los elementos. El cálculo de ${\textstyle \alpha }$ a partir de ${\textstyle F_{s}}$ es empírico y proviene del interés de dar una gama válida de valores de transición.

Figura 51: En la Figura (A), el elemento 4 tiene el mismo tamaño que el 1 y el 2 porque su nueva altura se calcula sólo a partir del $l$. En la Figura (B), $h_{2}$ es mayor que $h_{1}$ porque $h_{2}$ se calcula a partir del tamaño óptimo almacenado en el elemento 3.

Otro método que da mejores transiciones es guardar para cada lado del frente su tamaño óptimo de nuevo elemento y calcular el tamaño ${\textstyle h}$ en (99) a partir de este tamaño guardado en lugar de a partir del tamaño real del lado. La mejora del segundo método respecto al primero proviene de que si en los casos que al avanzar el frente se aprovecha puntos existentes, las dos magnitudes se hacen distintas. En estos casos, mantener el tamaño deseado provoca una mejora en el tamaño óptimo. La Figura 51, compara la opción de guardar ese tamaño óptimo con la de no guardarlo.

(¹) En este caso, el concepto tamaño óptimo se refiere a la altura del elemento.

4.14.2 Elección del punto para la creación del nuevo elemento

Una vez calculado el punto óptimo, es necesario encontrar todos los puntos existentes cercanos para escoger cual de ellos debe ser usado. La comparación de distancias de todos los puntos del frente es extremadamente costosa en tiempo y operaciones y daría lugar a algoritmos de orden ${\textstyle O^{2}}$ o peores. Es necesario, por tanto, el uso de algún algoritmo con una optimización mayor en la búsqueda de puntos alrededor de uno dado. El algoritmo usado es el que se basa en la estructura de almacenamiento de tipo Quadtree para dos dimensiones u Octtree para tres dimensiones. Estos algoritmos se describirán en el Apartado 4.12.1.

Como ya se comentó en el Apartado 4.12.1.3, para que los siguientes algoritmos funcionen correctamente, es más importante fijar un número aproximado de puntos de búsqueda, que definir un radio y su esfera asociada.

4.14.3 Borrado de elementos

En algunos casos, especialmente en la generación de malla de tetraedros, es posible llegar a una configuración del frente que no permita continuar en el avance. La razón de esta inestabilidad en el método proviene del hecho que la decisión de crear nuevos elementos proviene de la suma de un conjunto de tolerancias no directamente relacionadas entre ellas. Dada la construcción de un nuevo punto a una distancia determinada de una cara existente, aceptable según las tolerancias predefinidas, puede suceder que, más tarde, al proceder en el avance, la inserción de un nuevo elemento en el espacio dejado entre ese punto y la cara no cumpla con los requisitos establecidos. Hablando con más generalidad, puede producirse una acumulación de estado distorsionado mediante la adición de varias tolerancias mínimas que conlleve una situación de no solución. En la Figura 52, podemos ver dos tetraedros que no permiten la inserción de un nuevo elemento entre ellos.

Figura 52: Configuración sin posibilidad de continuar la generación de malla.

Cuando se llega a una situación de este tipo, se puede proceder de dos maneras, o bajando las exigencias de calidad de los elementos, posibilidad que se describirá en el Apartado 4.14.4, o borrando algunos elementos ya generados y reanudando la generación con un nuevo frente.

El proceso de borrado se realiza en dos fases:

1. Borrado corto.
2. Borrado de región.

El borrado corto consiste en borrar todos los elementos cuyas caras intersecten con las caras del elemento que se está intentando crear. En este caso se vuelve a intentar la generación desde la misma cara.

El borrado de región se realiza cuando ha fracasado el anterior y consiste en borrar todos los elementos asociados a las caras pertenecientes a la región de interés. En este caso, se realizan varios intentos de generación partiendo desde diferentes caras cuyo tamaño no difiera en exceso de la cara original.

4.14.4 Disminución de tolerancias

Otra técnica para solventar situaciones degeneradas consiste en la disminución de tolerancias. La idea es que en el avance normal, se requiere una calidad alta de elemento creado. Además también se requiere que la distancia a otras entidades, como a una cara ya existente en el frente, sea suficientemente grande.

En concreto, un nuevo punto debe cumplir que:

${\vec {N}}\cdot ({\vec {P}}-{\vec {C}}_{c})>\epsilon _{r}\cdot h$

donde ${\textstyle {\vec {N}}}$ es la normal al elemento, ${\textstyle {\vec {P}}}$ es el punto que se está comprobando, ${\textstyle {\vec {C}}_{c}}$ es el centro de la cara y ${\textstyle h}$ es una magnitud característica del nuevo elemento como su altura.

Por otra parte, las nuevas caras que se van a crear y que unen el nuevo punto con la cara del frente activo de la que se parte, deben cumplir una serie de comprobaciones de intersección con caras y aristas existentes. En la intersección cara-arista, se considera a la cara como un triángulo con funciones de forma ${\textstyle N(\alpha ,\beta )}$ y a la línea como una función del tipo ${\textstyle N(t)}$ con ${\textstyle t\in [0,1]}$. La intersección se produce si existe un punto en el espacio de parámetros ${\textstyle (\alpha ,\beta ,t)}$ tal que los dos puntos calculados como:

${\vec {P}}_{1}=(1-N(t)){\vec {P}}_{L1}+N(t){\vec {P}}_{L2}$ (101) ${\vec {P}}_{2}=\sum _{i=1}^{3}N_{i}{\vec {P}}_{Ti}$ (102)

coinciden. Esta sería la intersección exacta.

Para el cálculo de estas intersecciones en el interior del algoritmo de generación frontal, se usa la intersección extendida, que consiste en:

$(1-N(t)){\vec {P}}_{L1}+N(t){\vec {P}}_{L2}=\sum _{i=1}^{3}N_{i}(\alpha ,\beta ){\vec {P}}_{Ti}$ (103) se considera intersección extendida si: $\alpha ,\beta ,t\in [-\epsilon _{r},1+\epsilon _{r}]$ (104)

Considerando ${\textstyle \epsilon _{r0}}$ como la tolerancia adimensional definida, cuando se llega a una situación donde no se puede avanzar, se procede a disminuir este parámetro a ${\textstyle \epsilon _{r1}}$ y proseguir el avance. En caso de nuevo fallo, se intenta una vez más con un ${\textstyle \epsilon _{r2}}$ aún menor. En caso de fallar de nuevo se continúa con otros métodos como el de borrado de elementos definido en el apartado anterior. Si mediante una de estas reducciones de calidad se consigue generar un nuevo elemento, la generación prosigue con ${\textstyle \epsilon _{r0}}$ para que los elementos siguientes sean de alta calidad.

4.14.5 Ordenación de intentos

En el esquema de la Figura 53, se esquematiza el proceso de disminución de tolerancias y de borrado de elementos según el orden en que se van aplicando.

Figura 53: Orden de aplicación de los algoritmos de recuperación de errores.

4.15 MÉTODO DE DELAUNAY

El método de Delaunay es una técnica alternativa para la generación de malla no estructurada triangular o tetraédrica en un dominio. En [6] se describen detalladamente todos los algoritmos asociados con este método. Siguiendo el método constructivo que se describirá en siguientes apartados, se obtiene una malla de Delaunay que debe cumplir la siguiente definición:

Definición: Una malla triangular será malla Delaunay si cualquier círculo circunscrito a un triángulo no contiene en su interior ningún otro nodo de la malla¹.

La descripción del método que sigue está centrada en mallas triangulares. Para mallas volumétricas se podría crear una metodología similar.

Figura 54: Criterio de aceptación de un elemento según el esquema de Delaunay.

Dado un contorno plano, el proceso constructivo será el siguiente:

1. Se calcula una región rectangular que contenga al dominio.
2. Se generan nodos en el interior de esta región que cumplan los criterios de espaciamiento requeridos. A los nodos generados se les unen los nodos del contorno. Las técnicas de generación de puntos se describen en [6].
3. Se van comprobando nodos cercanos tres a tres. Se crea el círculo definido por tres nodos, si este círculo no contiene a ningún otro nodo se acepta el elemento. En la Figura 54, se visualiza este criterio. Los elementos que contienen puntos del contorno deben comprobarse especialmente pues es posible que el segmento de unión entre dos nodos del contorno no forme un triángulo que cumpla el criterio.
4. Cuando todo el rectángulo exterior está completamente cubierto de triángulos, es necesario recuperar los elementos que están en el interior del contorno. Para ello se empieza a buscar elementos a partir de un segmento de contorno cualquiera y se va avanzado mediante la búsqueda de elementos vecinos. Para ello debe mantenerse información de elementos vecinos, así como de los contornos.
5. Se aplican técnicas de alisado similares a las de los otros métodos.

(¹) Para malla de tetraedros puede usarse una definición similar sin más que substituir círculos por esferas.

4.15.1 Ventajas e inconvenientes del Método de Delaunay

Las ventajas e inconvenientes de unos métodos respecto a otros pertenecen en gran manera al campo subjetivo. De todas maneras, se puede intentar un pequeño resumen que de una idea a la persona interesada en el tema.

4.15.1.1 Ventajas

• Gran velocidad por el menor número de operaciones. El número de operaciones por elemento es menor debido a que en la creación del elemento por este método sólo se compara círculos con puntos, en contraposición al avance frontal donde se realiza el cálculo de gran número de intersecciones.
• Siempre se obtiene una malla. Independientemente de la calidad de los elementos, siempre se llega al final con este método, en contraposición al avance frontal donde puede llegarse a situaciones de no avance.
• Especialmente útil en el caso de que ya se disponga a priori de los nodos de la futura malla y únicamente quieran crearse los elementos.

4.15.1.2 Inconvenientes

• Control limitado en las transiciones de tamaños de elementos. Dificultad, en general, de crear unos puntos que cumplan con unos requerimientos de separación determinados. En el método de avance frontal, este requerimiento se impone de manera natural al crear los elementos.
• Mala calidad de los elementos resultantes. Esto se da especialmente en 3D para tetraedros. Es habitual obtener elementos de volumen nulo o negativo que es necesario corregir con el alisado.
• Dificultades adicionales en los puntos del contorno por la necesidad de mantener las caras del contorno, que son una información invariable de partida.

4.16 ALISADO

Se denomina alisado a las técnicas de mejora de calidad de la malla. Habitualmente estos algoritmos se aplican automáticamente sobre la malla recién generada. Puede considerarse, por tanto, como parte integrante del proceso de generación de malla.

Para entender lo que significa mejorar la calidad de la malla, hay que precisar que se entiende por calidad de malla. Genéricamente hablando, se considera la calidad de la malla en relación a la calidad de los elementos que la componen y calidad de los elementos, a la semejanza que tengan con un elemento regular, sea un triángulo equilátero o un tetraedro regular.

4.16.1 Indicador de calidad $I_{c}$4.16.1 Indicador de calidad I c {\displaystyle I {c}}

Es necesario definir un indicador de calidad que proporcione una medida de que elementos son mejores que otros y de cuales son aceptables y cuales no lo son. Existen multitud de indicadores descritos en la literatura (ver [9]). Algunos ejemplos de ellos son:

1. Relación de tamaños cara máxima/mínima-volumen.
2. Criterios de ángulos mínimos/máximos.

Se ha adoptado para todos los cálculos el criterio del coseno del ángulo mínimo. Dicho criterio se define de la siguiente manera:

Figura 55: Indicador de calidad $I_{c}$ para triángulos y para tetraedros.

Dados los puntos ${\textstyle {\vec {P}}_{1},{\vec {P}}_{3},{\vec {P}}_{3}}$ definidores del triángulo de estudio, se llama indicador de calidad ${\textstyle I_{c}}$ a:

$I_{c}=\min \left({\frac {({\vec {P}}_{i-1}-{\vec {P}}_{i})\cdot ({\vec {P}}_{i+1}-{\vec {P}}_{i})}{\|{\vec {P}}_{i-1}-{\vec {P}}_{i}\|\cdot \|{\vec {P}}_{i+1}-{\vec {P}}_{i}\|}}\right)\qquad \left|{\begin{array}{l}i=1,2,3\\i_{0}=i_{3}\\i_{4}=i_{1}\end{array}}\right.$ (107) si la normal a la superficie es ${\vec {N}}_{s}$, la normal al triángulo ${\vec {N}}_{t}$ es: ${\vec {N}}_{t}={\frac {({\vec {P}}_{i-1}-{\vec {P}}_{i})\times ({\vec {P}}_{i+1}-{\vec {P}}_{i})}{\|{\vec {P}}_{i-1}-{\vec {P}}_{i}\|\cdot \|{\vec {P}}_{i+1}-{\vec {P}}_{i}\|}}$ (108) entonces debe cumplirse que: ${\vec {N}}_{s}\cdot {\vec {N}}_{t}>0$ (109)

La condición (109) debe cumplirse para asegurar que el elemento no ha quedado mal orientado¹ En el caso de malla bidimensional, se supone al elemento situado en el plano ${\textstyle z=0}$ y la normal a la superficie será ${\textstyle {\vec {N}}=[{\begin{array}{ccc}0&0&1\end{array}}]^{T}}$.

Para tetraedros, el ángulo a considerar será el ángulo diedro entre los planos de dos caras. Si llamamos ${\textstyle {\hat {N}}_{i}}$ a las normales unitarias a las caras, el indicador ${\textstyle I_{c}}$ se calculará como:

$I_{c}=\min({\hat {N}}_{i}\cdot {\hat {N}}_{j})\qquad \left|{\begin{array}{l}i\neq j\\i,j=1,\dots ,4\end{array}}\right.$

(110)

Para asegurarnos que el volumen del elemento no resulta negativo, puede calcularse su Jacobiano, que debe ser positivo. En general no es necesario este cálculo, debido al hecho de que el método de generación impedirá la creación de elementos degenerados.

Se observa que este indicador de calidad será mayor como peor sea la calidad y en su extremo, cuando el ángulo sea nulo, el indicador valdrá 1.

(¹) A priori, el método de generación debe impedir la creación de elementos mal orientados o de volumen negativos pero pueden darse casos degenerados.

4.16.2 Proceso de mejora de la malla

Dada una malla a la que se desea hacer un proceso de mejora, sería posible aplicar los algoritmos que se van a describir en siguientes apartados a todos los elementos y nodos que la forman. Este criterio daría mallas de calidad alta pero sería absolutamente prohibitivo a nivel de coste computacional en el momento que las mallas se hiciesen un poco grandes. Por otro lado, interesa iterar varias veces sobre los algoritmos porque las mejoras de calidad en los elementos no son acumulativas, sino que tienen interdependencia unas con otras.

Asimismo, la pretensión de estos procesos no es tanto la mejora global de la malla como la mejora de los elementos peores. Éstos se consideran como críticos para el éxito de cualquier tipo de análisis posterior con esta malla.

El proceso será calcular el índice ${\textstyle I_{c}}$ para todos los elementos de la malla y escoger los peores en un porcentaje prefijado. Seguidamente se aplicará uno de los algoritmos de mejora y se volverá a rehacer la lista de los peores. Se aplicará un nuevo algoritmo y se procederá de igual manera. Es conveniente iterar dos o tres veces sobre el conjunto de los algoritmos para obtener una buena relación mejora de calidad-tiempo de computación.

4.16.3 Algoritmos de mejora de la malla

Se van a presentar dos algoritmos de mejora de la calidad de la malla desarrollados en la tesis. Uno de ellos se puede aplicar únicamente a triángulos y el otro tanto a triángulos como a tetraedros. Los algoritmos son:

• Intercambio de diagonales.
• Movimiento de nodos en el macro-elemento asociado.

Ambos algoritmos parten de un elemento, que es el que se considera que se debe mejorar.

4.16.3.1 Intercambio de diagonales

Este algoritmo es únicamente válido para triángulos. Dado el elemento a mejorar, se va comparando alternativamente con las caras que tiene a su alrededor, que serán 1,2 ó 3 caras, y se basa en que, dados los dos triángulos que comparten una cara, se comprueba sus indicadores de calidad ${\textstyle I_{c1}}$ y ${\textstyle I_{c2}}$. Si los triángulos son ${\textstyle {\vec {P}}_{1},{\vec {P}}_{2},{\vec {P}}_{3}}$ y ${\textstyle {\vec {P}}_{2},{\vec {P}}_{4},{\vec {P}}_{3}}$ según la Figura 56, se substituyen dichos triángulos por los ${\textstyle {\vec {P}}_{1},{\vec {P}}_{2},{\vec {P}}_{4}}$ y ${\textstyle {\vec {P}}_{1},{\vec {P}}_{4},{\vec {P}}_{3}}$. Se calcula de nuevo los indicadores de calidad y si ambos son mayores que los anteriores, se acepta el cambio.

Figura 56: Algoritmo de intercambio de diagonales.

Debe tenerse en cuenta que no interesa que el conjunto de ambos elementos mejore, sino que sólo aceptaremos el cambio si los dos nuevos elementos mejoran a los dos antiguos. La razón de ello estriba en el hecho de que el objetivo es mejorar los elementos peores de la malla en lugar de querer mejorar la calidad media.

La razón de la no aplicabilidad de este algoritmo a tetraedros viene dada por el hecho de que si se substituye una cara por la cara asociada, se cambia la conectividad con otros tetraedros vecinos y se pierda la conformidad de la malla. Esta situación se ve claramente en la Figura 57.

Figura 57: Si cambiamos la cara $\bigtriangleup _{235}$ por la cara $\bigtriangleup _{134}$, los dos nuevos tetraedros son válidos, pero al hacerlo cambiamos también la conectividad de las caras $\bigtriangleup _{125}$ y $\bigtriangleup _{245}$ hacia $\bigtriangleup _{124}$ y $\bigtriangleup _{145}$, lo cual dará no conformidad con los elementos vecinos.

4.16.3.2 Movimiento de nodos en el macro-elemento asociado

Otro criterio de mejora de la malla consiste en que, dado el elemento a mejorar, se aplica el siguiente algoritmo a cada uno de sus nodos que no estén en el contorno:

Se buscan todos los elementos que contienen a ese nodo. Se calcula el índice de calidad ${\textstyle I_{c}}$ para todos los elementos y nos quedamos con el mínimo ${\textstyle I_{c,min0}}$. Se realiza un movimiento del nodo según alguna de las técnicas descritas posteriormente y se calcula de nuevo ${\textstyle I_{c,min1}}$ para el conjunto de los elementos. Si ${\textstyle I_{c,min1}<I_{c,min0}}$ se acepta el cambio y se procede hacia un nuevo intento de movimiento.

Puede observarse que, al igual que en el algoritmo de intercambio de diagonales, es preferible mantener el ${\textstyle I_{c}}$ del mínimo que aumentar la media de calidad del conjunto.

Las dos técnicas de movimiento de nodos desarrolladas son:

1. Movimiento hacia el centroide del conjunto de elementos.
2. Movimiento de factor según las diagonales del conjunto.

Dichas técnicas se describen a continuación.

Figura 58: En la Figura (A), el movimiento de nodo se produce hacia el centroide del conjunto. En la (B), se produce un movimiento de factor en una de las diagonales del conjunto.

4.16.3.3 Movimiento hacia el centroide del conjunto de elementos

El movimiento queda definido por:

${\vec {P}}_{cn}={\frac {\sum _{i=1}{n}{\vec {P}}_{i}}{n}}$

donde ${\textstyle {\vec {P}}_{cn}}$ es la nueva posición del nodo y ${\textstyle {\vec {P}}_{i}}$ son los ${\textstyle n}$ nodos que rodean al nodo de estudio. En la Figura 58 (A), se representa gráficamente esta translación del nodo.

Este esquema se utiliza principalmente en mallas de triángulos, pues en la mayoría de los casos mejora la calidad del conjunto a un nivel muy cercano al óptimo. En cambio, en el caso de tetraedros, es más frecuente que esta nueva posición origine elementos de peor calidad y que incluso, la posición final del nodo quede fuera del conjunto. En este último caso quedarían elementos de área negativa.

4.16.3.4 Movimiento de factor según las diagonales del conjunto

La técnica consiste en mover el nodo de estudio en una pequeña magnitud tal que el punto final se mantenga en el interior del conjunto de elementos y que mejore la calidad del elemento peor.

Una manera de definir las direcciones de movimiento viene dada por escoger las diagonales que unen el nodo con todos los otros nodos que le rodean. Es una elección arbitraria pero cómoda, ya que obtenemos con poca cantidad de cálculos un conjunto de direcciones en el espacio que tienden a repartirse homogéneamente.

La magnitud del movimiento ${\textstyle \alpha L}$ se define de manera relativa a la longitud de la diagonal ${\textstyle L}$, y se va probando de manera iterativa para diferentes valores predefinidos del factor ${\textstyle \alpha }$. La posición final vendrá dada por:

${\vec {P}}_{cn}=\alpha {\vec {P}}_{i}+(1-\alpha ){\vec {P}}_{cv}$

siendo ${\textstyle {\vec {P}}_{i}}$ uno de los nodos que rodean al nodo de estudio.

4.17 GENERACIÓN DE MALLA DE SUPERFICIE EN EL ESPACIO

En el Apartado 4.11, se describían las técnicas de generación que requerían una transformación previa al plano. Ahora se expondrán sus equivalentes en caso de generación en el espacio.

Generación en el espacio significa que la implementación de cualquier método se realiza directamente en la posición que ocupa la geometría original. El método usado para este tipo de generación es el de generación frontal. Las diferencias con el método tradicional se describirán a continuación.

Los pasos y particularidades de la generación frontal en su modalidad espacial son:

1. Se elige la cara del frente donde se va a continuar el avance, de igual forma que en la generación 2D.
2. El punto óptimo del nuevo triángulo a generar se calcula en el plano local tangente. Este es el plano ortogonal a la normal a la superficie en la arista de avance.
3. El punto debe resituarse sobre la superficie desde el plano local. Para ello se realiza una nueva proyección sobre la superficie.
4. La comprobación de intersecciones debe realizarse mediante la comparación de caras en el espacio.
5. La comprobación de la calidad y aceptabilidad de un elemento debe hacerse siempre en relación a la normal de la superficie en ese punto.
6. Al realizar el movimiento de nodos en el proceso de mejora de la malla debe siempre resituarse los nodos sobre la superficie original.

A continuación se describirán estas etapas con mayor detalle.

4.17.1 Cálculo del punto óptimo

El cálculo del punto óptimo para el nuevo elemento que se va a generar se realiza de manera equivalente al caso bidimensional pero situado encima del plano tangente. Para ello es necesario conocer previamente la normal a la superficie. Como se conocen los parámetros ${\textstyle u,v}$ para todos los nodos existentes, el problema se reduce al cálculo de la normal a la superficie en ${\textstyle u=(u_{1}+u_{2})/2,v=(v_{1}+v_{2})/2}$ siendo ${\textstyle u_{1},v_{1}}$ las coordenadas paramétricas del nodo 1 de la arista y ${\textstyle u_{2},v_{2}}$ las del nodo 2.

El punto calculado no residirá en el interior de la superficie. Es necesario hacer una operación de reposicionamiento sobre ella (ver el Apartado 3.10.3). En la Figura 59, se puede ver una descripción gráfica del proceso.

Figura 59: Cálculo del punto óptimo del nuevo elemento a generar.

4.17.2 Intersecciones entre caras

Al referirnos a generación de malla espacial no es posible hablar propiamente de intersección de caras, pues éstas no intersectarán directamente en multitud de casos. Más bien, debemos referirnos al concepto de intersección de caras proyectadas. En la Figura 60, se aprecia un caso de este tipo.

Para considerar que dos caras en el espacio intersectan, debemos referirnos de nuevo al plano tangente. Se realiza una transformación de los nodos que forman las caras a un nuevo sistema de coordenadas tal que su eje ${\textstyle z}$ coincide con la normal a la superficie.

${\vec {P}}_{i,L}=T^{t}P_{i}$ $\Delta z=\|I_{1,z,L}-I_{2,z,L}\|$

(111)

El cálculo de la intersección se realiza sin tener en cuenta la coordenada ${\textstyle z}$ local de los triángulos. Los puntos de intersección línea-triángulo ${\textstyle I_{i}}$ tendrán unas coordenadas ${\textstyle z}$ diferentes para cada uno de los triángulos. La magnitud ${\textstyle \Delta z}$ será la diferencia entre coordenadas ${\textstyle z}$ locales de esos puntos. Se considerará intersección si existe algún ${\textstyle I_{i}}$ tal que ${\textstyle \Delta z}$ sea menor que una magnitud relativa dada.

La razón de tener que limitar las distancias ortogonales al plano de proyección en las intersecciones, viene dada para distinguir casos anómalos de altas curvaturas en la superficie.

Figura 60: Diversas vistas de dos caras que no intersectan en el espacio pero que deben considerarse como una intersección para el algoritmo de avance frontal.

4.17.3 Control de calidad de los elementos

En el plano existe una orientación dada de manera implícita que consiste en aplicar la regla dextrógira a una normal definida en sentido ${\textstyle z}$ positiva. Ello resulta en un sentido positivo de los giros antihorarios.

En el espacio es necesaria una referencia adicional a la orientación de la superficie. Al indicador ${\textstyle I_{c}}$ definido en el Apartado [[#4.16.1 Indicador de calidad $I_{c}$|4.16.1]], hay que añadir un criterio adicional de comparación de normales. En concreto, para un triángulo ${\textstyle {\vec {P}}_{1},{\vec {P}}_{2},{\vec {P}}_{3}}$ se puede expresar:

${\vec {N}}_{t}=({\vec {P}}_{2}-{\vec {P}}_{1})\times ({\vec {P}}_{3}-{\vec {P}}_{1})$ $I_{a}={\vec {N}}_{t}\cdot {\vec {N}}_{s}$

(112)

con ${\textstyle {\vec {N}}_{s}}$ la normal a la superficie. Si ${\textstyle I_{a}}$ es menor que 0 debe considerarse el elemento mal orientado y por tanto inaceptable.

4.17.4 Comparación de la generación en el plano versus en el espacio

La operación básica que distingue la generación en el plano en relación a la generación en el espacio es la de reposicionamiento de un punto sobre la superficie.

Esta operación es altamente costosa en tiempo de computación y puede llegar a no obtener un resultado correcto para superficies muy distorsionadas. Aunque también es necesario realizarla en la generación de malla plana para los contornos, en la generación de malla espacial se usa mucho más extensamente. Por otro lado, esta operación puede especializarse en los contornos de tal manera que tenga resultados más precisos.

Figura 61: Velocidad en elementos por minuto con el algoritmo de generación en el plano (pentium 200). El gráfico total representa el tiempo total de generación. El AF representa el tiempo del avance frontal.

Existen otros condicionantes en la generación espacial que dificultan la obtención de una malla. Entre ellos se puede mencionar la dificultad de atravesar zonas de alta curvatura con elementos grandes o la confusión que provocan los contornos cerrados.

El resultado de ello es que la generación de malla plana es mucho más rápida que la homóloga en el espacio y es más fiable. Como contrapartida, puede dar elementos más distorsionados que la generación espacial.

Figura 62: Velocidad en elementos por minuto con el algoritmo de generación en el espacio (pentium 200). El gráfico total representa el tiempo total de generación. El AF representa el tiempo del avance frontal.

En las Figuras 61 y 62, puede observarse que la velocidad de generación plana es un orden de magnitud mayor que la espacial. En todo el proceso de generación hay algunos algoritmos que son de orden lineal y otros que son de orden ${\textstyle n\log(n)}$. A partir de las gráficas se deduce que después de un incremento de la velocidad al inicio debido a tareas de coste fijo, se entra en una curva ligeramente descendente pero bastante constante para el rango de elementos analizado. Se deduce de ello que los algoritmos de mayor coste computacional son los de orden lineal.

Figura 63: Mallas obtenidas en una superficie distorsionada según la generación plana y la espacial.

Figura 64: Gráfica de calidad de elementos para una superficie distorsionada según la generación plana y la espacial.

4.17.5 Tiempos relativos de computación

Number of Elems=51775

Time Table: ...............  cputime[s] (percent[%]) Uses – Elems/min

.............................. 18.72000[s] (100.00000[%]) 1 – 165946
. beginadvancingfront ........... 12.87000[s] (68.75000[%]) 1 – 241375
.  . Advancing front ............... 11.17000[s] (86.79099[%]) 1 – 278111
.  . Smoothing ..................... 1.44000[s] (11.18881[%]) 1 – 2157292
.  .  . Calc Quality .................. 1.03000[s] (71.52778[%]) 4 – 3016019
.  .  . Swap Diagonals ................ 0.06000[s] (4.16667[%]) 1 – 51775000
.  .  . Move Points To Center ......... 0.35000[s] (24.30556[%]) 2 – 8875714
.  . Other operations .............. 0.26000[s] (2.02020[%]) – 11948077
. Other operations .............. 5.85000[s] (31.25000[%]) – 531026

En las Tablas 4 y 5, se puede observar el porcentaje de tiempo que se dedica a cada una de las tareas en la generación de malla por el método del avance frontal. A partir de ellas se pueden deducir algunas conclusiones:

• El porcentaje de tiempo dedicado a la mejora de la malla o alisado es del 10 al 14% Este número se consigue mantener en esa magnitud debido a que la mejora sólo se realiza en los elementos de peor calidad. Debido a este criterio, la mayoría de elementos no son modificados. Así conseguimos el objetivo de mejorar la calidad mínima de la malla en lugar de mejorar la calidad media.
• Dentro de los algoritmos de mejora, el que más tiempo consume es el de movimiento de nodos.
• Los algoritmos de borrado de elementos en la generación de tetraedros ocupan, para este ejemplo, un tiempo despreciable. En otros modelos más complejos, este tiempo puede aumentar pero siempre será pequeño en relación al tiempo total.
• La mayor parte del tiempo de computación para malla tridimensional, se usa en encontrar los nodos cercanos a la cara activa, a partir de ellos obtener las aristas y caras cercanas, y finalmente, escoger el elemento óptimo mediante la comprobación de intersecciones y otros criterios.

Number of Elems=14854

Time Table: ................. cputime[s] (percent[%]) Uses – Elems/min

Absolute total................ 49.02[s] (100.00000[%]) 1 – 18181
+ Advancing front total ......... 47.21[s] ( 96.30763[%]) 7 – 18878
|  + Advancing front ............. 40.36[s] ( 85.49036[%]) 7 – 220826
|  |  + CalcBestPosition ............  0.10[s] (  0.24777[%]) 14857 – 8912400
|  |  + FindClosePoints3D ...........  9.86[s] ( 24.43013[%]) 14857 – 90389
|  |  + CalcCloseFaces3D ............  5.42[s] ( 13.42914[%]) 14857 – 164435
|  |  + CalcCloseAristes3D .......... 10.12[s] ( 25.07433[%]) 14857 – 88067
|  |  + FindBest3D .................. 12.28[s] ( 30.42616[%]) 15182 – 72576
|  |  + EnterFace ...................  0.19[s] (  0.47076[%]) 28835 – 4690736
|  |  + TakeFaceOut .................  0.20[s] (  0.49554[%]) 30589 – 4456200
|  |  + DeleteClosedElemsShort ......  0.02[s] (  0.04955[%]) 1 – 44562000
|  |  |  + EnterFace ................... 0.00[s] (  0.00000[%]) 3
|  |  |  + TakeFaceOut ................. 0.00[s] (  0.00000[%]) 5
|  |  + Other operations ............ 2.17[s] (  5.37661[%])   – 410709
|  + Smoothing ................... 6.77[s] ( 14.34018[%]) 7 – 131645
|  |  + Calc Quality ................ 0.50[s] (  7.38552[%]) 27 – 1782480
|  |  + Swap Diagonals .............. 0.02[s] (  0.29542[%]) 6 – 44562000
|  |  + Move Points To Center ....... 6.25[s] ( 92.31905[%]) 14 – 142598
|  + EnterFace ................... 0.01[s] (  0.02118[%]) 1756 – 89124000
+ Other operations ............ 1.81[s] (  3.69237[%])   – 492397

4.18 GENERACIÓN DE MALLA EN LAS SUPERFICIE-MALLA

En el Apartado 3.10.7, llamamos Superficie-malla a un nuevo tipo de entidad geométrica definida mediante una malla de superficie. Para poder generar sobre dicha superficie, se usa el algoritmo de generación de malla en el espacio con algunas variaciones. Los datos que la entidad geométrica debía proporcionar al algoritmo de generación eran:

• Dado un punto en la superficie, calcular la normal a la superficie en ese punto.
• Dado un punto en el espacio, en una posición cercana a la superficie, encontrar el punto perteneciente a dicha entidad que cumpla el criterio de mínima distancia. O sea, hacer el reposicionamiento de dicho punto.

El proceso general para calcular estos datos consiste en encontrar al elemento en la malla de soporte al cual pertenece el punto¹.

(¹) Un punto pertenecerá a un elemento cuando esté en su interior o cuando su proyección esté en el interior de dicho elemento.

4.18.1 Búsqueda del elemento activo

La búsqueda del elemento activo de malla soporte se realiza mediante el criterio del camino (ver [20]). Según este criterio y dado un elemento cualquiera, es posible calcular las coordenadas naturales más la altura de ese punto respecto al elemento.

El proceso es el siguiente: Dado un triángulo o cuadrilátero, se calcula su normal ${\textstyle {\vec {N}}}$ y se define una matriz de transformación ${\textstyle T}$ tal que el eje ${\textstyle z'}$ de los ejes locales que definen dicha rotación coincida con la normal ${\textstyle N}$. Los nodos que definen el triángulo y el punto de referencia pueden convertirse a ese sistema local mediante:

${\vec {P}}_{i}'=T^{t}\cdot {\vec {P}}_{i}$

En ese sistema la ${\textstyle z}$ del punto será ${\textstyle z=\|{\vec {P}}_{r,z}-P_{1,z}\|}$ las coordenadas naturales ${\textstyle L_{1},L_{2}}$ y ${\textstyle L_{3}}$ se calculan resolviendo el sistema:

$P_{r,x}$ $=L_{1}P_{1,x}+L_{2}P_{2,x}+L_{3}P_{3,x}$ $P_{r,y}$ $=L_{1}P_{1,y}+L_{2}P_{2,y}+L_{3}P_{3,y}$ $L_{1}$ $=1-L_{2}-L_{3}$

Si ${\textstyle L_{1},L_{2},L_{3}\in [0,1]}$, el punto pertenece al triángulo. La ${\textstyle z}$ del punto debe limitarse para distinguir el caso de un punto que proyecte en el interior de dos triángulos que pertenecen a dos zonas distintas de la malla o que tienen el ángulo entre ambos es muy grande.

Figura 65: Búsqueda del nuevo elemento en el algoritmo del camino en función de las coordenadas naturales.

Para encontrar el elemento, se parte de un elemento cualquiera. Habitualmente se usa el que se ha obtenido en un paso anterior pues la probabilidad de que sea cercano al que buscamos es alta. Se calculan las coordenadas naturales locales del punto en relación a este elemento. Si alguna de ellas está fuera del intervalo ${\textstyle [0,1]}$, nos indicará hacia que lado del triángulo se debe avanzar. En la Figura 65, se puede observar el criterio de lado preferido. Como se ha almacenado información de elementos vecinos, se puede realizar el mismo proceso para el nuevo triángulo y proceder, de esta manera, en forma iterativa. Se debe mantener una lista de elementos ya visitados y fijar un número máximo de elementos en el camino para no entrar en recursiones infinitas.

En algunos casos, no es posible llegar al elemento buscado mediante este procedimiento. Entonces, se recurre al método de buscar la arista más cercana al punto y escoger, de entre los dos elementos que contienen a la arista, al que tenga unas coordenadas naturales mejores.

4.18.2 Cálculo de la nueva normal

Una vez obtenido el elemento óptimo, la normal se puede calcular como la normal al elemento. Para conseguir geometrías más suaves, puede calcularse un campo de normales mejorado. Este pasa por calcular previamente las normales en los nodos de la malla y, a partir de ellos, recalcular las nuevas normales en los elementos. De esta manera, la normal a un elemento tiene información de los elementos vecinos y se consigue más continuidad entre ellos.

Figura 66: Elección de la normal correcta para puntos situados en aristas.

Debe tenerse en cuenta que para un punto situado sobre una arista, la normal según uno u otro elemento de ambos lados de la arista serán muy diferentes si los elementos tienen un alto ángulo entre ellos. Esto será especialmente importante para el caso de aristas pertenecientes a líneas interiores del dominio (ver el Apartado 3.10.7). La elección del elemento correcto para este cálculo, requiere de información adicional. Esta información puede ser un segundo punto que nos de la orientación de avance en la generación de nuevos elementos. Una situación de este tipo puede verse en la Figura 66.

4.18.3 Reposicionamiento del punto

Cuando se ha obtenido el elemento óptimo y sus coordenadas naturales locales mediante la técnica descrita anteriormente, se considera al punto reposicionado como la transformación al espacio global del punto de coordenadas naturales equivalentes.

${\vec {P}}_{n}'$ $=\sum _{i=1}^{3}L_{i}{\vec {P}}_{i}$ ${\vec {P}}_{n}$ $=T\cdot {\vec {P}}_{n}'$

4.19 APORTACIONES EN LA GENERACIÓN DE MALLA

En este apartado se pretende distinguir los algoritmos y métodos que se han ido desarrollando a lo largo de este trabajo en contraposición al estado del arte y a los conocimientos acumulados sobre el tema. Más que pretender afirmar lo novedoso de algunas técnicas, se va a intentar mostrar cómo la utilización de algoritmos conocidos en maneras un tanto distintas, puede llevar a unos mejor cumplimiento de los objetivos, o sea, la generación eficiente de malla de calidad para modelos geométricos arbitrariamente complejos.

En las siguiente líneas se van a repasar algunos de los conceptos descritos a lo largo del capítulo y se van a comentar desde este punto de vista.

4.19.1 Transiciones de tamaños

El proceso anterior a la generación de malla que consiste en corregir los tamaños óptimos de elementos para compatibilizarlos entre ellos fue descrito en el Apartado 4.7.2. La innovación de más interés ha consistido en aplicar esta transición de tamaño a mezclas de mallas estructuradas y no estructuradas para minimizar el número de elementos que describen correctamente una geometría. Debe tenerse en cuenta que este preproceso que se produce anteriormente a la generación es crucial para obtener éxitos en la generación en modelos geométricos arbitrariamente complejos. También es básico para incrementar la automatización y minimizar, por tanto, el trabajo humano para obtener una malla correcta.

4.19.2 Superficies NURBS estructuradas

Tradicionalmente, los programas de generación de malla estructurada permiten la generación en superficies de 3 ó 4 lados ya que estos son los que forman la estructura óptima para este proceso. En los desarrollos de esta tesis se ha incluido la posibilidad de que la superficie tenga más de 4 lados siempre que éstos se puedan describir topológicamente mediante una forma cuadrilátera en el espacio. El proceso consiste en hacer abstracción de las líneas externas y convertirlas interiormente en 4 líneas para proceder a la generación. Este proceso tiene la gran ventaja que permite que los tamaños en lados no se transfieran a través de todas las superficies del modelo, pudiendo de esta manera poner más o menos elementos según zonas. El inconveniente es que ha sido necesaria la creación de algoritmos de compatibilización de tamaños entre superficies que pueden llevar, en un caso extremo, a imposibilitar una malla conforme a lo largo de todo el dominio.

4.19.3 Generación de malla estructurada

La generación de malla estructurada se realiza, tanto en las superficies como en los volúmenes, mediante la aplicación del algoritmo de interpolación Coon parcialmente bicúbico. Aunque este algoritmo ha sido y es utilizado ampliamente en muchas aplicaciones (ver [1]), su uso en la generación de malla estructurada representa un nuevo acercamiento al problema de insertar nodos en el interior de un dominio.

La concentración de elementos para dichas mallas se realiza con un incremento según una progresión geométrica que se asigna y se calcula línea por línea (ver Apartado 4.8.5). De esta manera se puede tener una combinatoria mayor para colocar los elementos en forma óptima.

4.19.4 Corrección NURBS para líneas

Uno de los problemas más graves que aparecen en la generación de malla sobre NURBS, consiste en que éstas pueden tener parametrizaciones muy alejadas a la del parámetro arco. En generación de malla, es completamente imprescindible tener un buen control en los tamaños de elementos en cada zona. Para solventar este problema se ha procedido de dos maneras:

1. Corrección y modificación de las entidades geométricas como paso previo a la generación de malla.
2. Corrección iterativa en el cálculo del incremento del parámetro en el momento de calcular los nodos sobre las líneas NURBS.

La primera opción se describe en el capítulo sobre modelización geométrica. La segunda opción, descrita en el Apartado 4.9.2, consiste en una mejora importante para poder obtener una separación óptima entre los nodos de una línea. Ambas opciones son aportaciones de la tesis a la generación de malla sobre NURBS.

4.19.5 Concentraciones de elementos en líneas

Cuando hay que generar nodos sobre una línea con diferentes tamaños en diferentes zonas de dicha entidad, es necesario el uso de un algoritmo que posicione estos nodos de tal manera que se distribuyan correctamente a lo largo de todo ella. El algoritmo de la doble parábola con reposicionamiento iterativo desarrollado en la tesis cumple muy adecuadamente con dichos preceptos. En concreto, se adapta muy bien al parámetro de transición definido por el usuario.

4.19.6 Transformación al espacio equivalente en generación de malla de superficie

Para la generación de malla de superficie, las técnicas más habituales en la literatura consisten en la generación en el espacio paramétrico con corrección mediante la métrica del espacio. Algunos ejemplos de ello pueden verse en [8]. También se encuentran ejemplos en la literatura sobre generación en un espacio equivalente al paramétrico.

La definición de un espacio equivalente mediante las funciones de forma del cuadrilátero serendípito de 8 nodos desarrollada en la tesis, constituye un nuevo enfoque para esta técnica. Debe tenerse en cuenta que el cálculo de este cuadrilátero debe poderse realizar sobre superficies extremadamente distorsionadas, lo cual implica una alta sofisticación en la algoritmia.

4.19.7 Uso de Octtrees modificados

El uso de Octtrees para búsquedas optimizadas se ha convertido en habitual para todo generador de malla que pretenda ser eficiente computacionalmente [24]. Ha sido necesario en la tesis la implementación de funcionalidades adicionales para poderlo adaptar a todas las necesidades que se han ido presentado. En concreto:

• Se le ha dado capacidad de incrementar su tamaño dinámicamente respecto al Octtree inicial para poder tratar los casos en que no se conoce a priori el tamaño total del dominio.
• Cuando entran en la estructura un número de puntos coincidentes mayor que el número máximo de puntos que caben en un octante, se produce desbordamiento a no ser que se tomen medidas especiales para evitarlo. Las medidas desarrolladas consisten en que al disminuir el tamaño de un octante por debajo de una tolerancia dada, se cambia el proceso de inserción, convirtiendo a uno nuevo de orden cuadrático.
• El criterio de extracción de puntos cercanos cumpliendo los requisitos de obtención de una cantidad aproximadamente fija y de eficiencia computacional, han propiciado el desarrollo de un nuevo algoritmo complejo pero rápido para conseguir esos puntos óptimos.

4.19.8 Octtree unidimensional

La extensión, o en este caso más bien reducción, de la técnica Octtree al caso unidimensional, se ha usado para la ordenación eficiente de caras del frente de generación de malla.

Para conseguir una extracción rápida de las caras usadas, ha sido necesario desarrollar criterios adicionales de información cruzada. Debido al hecho de que muchas caras tienen tamaño equivalente, ha sido necesario implementar técnicas que no redujesen la ordenación a orden ${\textstyle O^{2}}$.

4.19.9 Generación frontal tridimensional

Los algoritmos de generación de tetraedros en el espacio han sido materia de investigación abundante en los últimos años. Ejemplos de literatura sobre el tema pueden encontrarse en [8] y [23].

Para aumentar la eficiencia del método, la calidad de los elementos resultantes y sobre todo, la seguridad en el éxito de la generación, se han añadido algunas técnicas a las ya descritas. En concreto:

• Se mantiene una historia del incremento del tamaño de los elementos, pues basar el incremento del tamaño óptimo únicamente en el tamaño actual del elemento del frente, evita que se produzca en muchos casos el crecimiento deseado. La solución consiste en que cada elemento contenga información no sólo de su tamaño, sino también del tamaño que debería haber tenido.
• El borrado de elementos, que ya ha sido descrito en la literatura (ver [3]), es una técnica eficiente cuando se llega a situaciones que impiden avanzar en el frente. En apartados anteriores se ha descrito la técnica desarrollada que consiste en combinar el borrado según zona de influencia de la cara activa con el borrado de una sub-zona mayor y intento de avance según diversas caras.
• El criterio de diferentes niveles de tolerancia desarrollado permite asegurar una alta calidad de la mayoría de los elementos y al mismo tiempo, permite proceder con una tolerancia menor en la generación cuando no es posible continuar con la óptima.
• La combinación de los dos algoritmos anteriores ha obligado a definir un esquema en la ordenación de intentos para poder proseguir en el caso de quedar un frente imposibilitado.

4.19.10 Mejora de la malla o alisado

Las técnicas de mejora de la malla han sido descritas ampliamente en la literatura de generación de malla (ver [7]). Más concretamente, tanto el algoritmo de intercambio de diagonales como el de movimiento de nodos en el centroide de superficies elementos se han convertido en estándares en la materia. En el algoritmo desarrollado en la tesis, la consecución de una buena mejora de la malla radica en la utilización de un indicador de error adecuado y en un buen criterio para reposicionar los nodos que se desplazan en el centroide.

4.19.11 Generación de malla de superficie en el espacio

La generación directamente en el espacio, con las dificultades del reposicionamiento de nodos sobre la superficie y las del cálculo de las pseudo-intersecciones de caras planas en el espacio ha permitido desarrollar un nuevo algoritmo que, aún siendo más lento que otros algoritmos alternativos, proporciona elementos de gran calidad para superficies altamente distorsionadas.

Por el hecho de que no puede llegarse siempre a una generación óptima, es imprescindible usar dicho algoritmo en combinación con otras técnicas de generación de malla en el plano. Aceptando todas estas dificultades, este método de generación es de gran utilidad en los casos en que la calidad de los elementos resultantes sea un requisito imprescindible para el método de cálculo.

4.19.12 Las superficies-malla

El método de generación de malla sobre un dominio superficial descrito mediante otra malla, ha sido ya implementado y descrito en [20]. Sin embargo, las dificultades de implementación práctica y eficiente requieren de todo un conjunto de técnicas que lo complementen.

En concreto, los puntos que provocan más dificultades a esta tecnología son:

• El salto en el frente a través de una línea de unión de elementos que defina una elevado ángulo entre las dos partes.
• La búsqueda eficiente del elemento que debe servir de base para el cálculo de normales o para el reposicionamiento del nuevo nodo requiere de métodos de avance según un camino a través de la malla y de cálculos rápidos de funciones de forma en coordenadas locales.

La aportación en las superficies-malla, ha consistido en la implementación eficiente de los algoritmos descritos.

4.20 CONCLUSIONES A LA GENERACIÓN DE MALLA

Durante los últimos años se ha producido una importante investigación, con su desarrollo asociado, de multitud de técnicas de generación de malla. Ello ha producido abundante literatura sobre el tema y la aceptación de algunas técnicas como las que deben ser usadas de facto, o sea, que han pasado ya al conocimiento general.

Sin embargo, la obtención a partir de unas técnicas genéricas de unos generadores de malla de robustez industrial y que sean lo suficientemente automáticos para casi anular la intervención del usuario, es aún motivo de importantes desarrollos y preocupaciones. Puede resumirse diciendo que es un tema muy avanzado a nivel de investigación pero que aún no ha llegado a su absoluta madurez en su aplicación práctica y eficiente.

Se espera que los desarrollos de esta tesis descritos en este capítulo, ayuden a avanzar la tecnología de generación de malla hacia ese objetivo de alta robustez y automatización.

5 Adaptación a un tipo de análisis

5.1 INTRODUCCIÓN

Se entiende por adaptación a un tipo de análisis al conjunto de todos los pasos, procesos y programación adicional que se requiere para que un sistema de preproceso pueda tratar y más tarde transferir toda la información necesaria para describir dicho análisis a un programa de cálculo.

Tal como se describía en el Capítulo 2, el sistema desarrollado en la tesis se ha organizado para poderse utilizar en cualquier tipo de análisis. Para poder cumplir tal objetivo, cabría pensar en tres posibilidades:

1. En la creación del sistema se incluyen todos los tipos de datos que se necesitarán para los análisis contemplados y también se incluye la adaptación para todos esos análisis.
2. En el interior del sistema se incluyen todos los tipos de datos que se necesitarán. Para cada tipo de análisis se escoge de entre esos datos.
3. Se considera a la adaptación como un proceso externo al sistema en sí que se le añade a posteriori. En cada adaptación con cada tipo de análisis se definen totalmente los datos a usar.

La primera posibilidad es la más cerrada en el sentido de que no admite ningún tipo de análisis nuevo. La segunda posibilidad permitiría incluir tipos de análisis nuevos siempre que no necesiten estructuras de datos nuevas. Así, si no se han contemplado los tipos de datos particulares que requiere un análisis concreto, la adaptación posterior será imposible.

La tercera posibilidad, que es la más genérica, es la que ha sido adoptada. Para hacerla posible, es necesario realizar una abstracción de todos los datos que se puedan necesitar para cualquier análisis y entenderlos como unas estructuras de datos más sencillas.

Con este criterio se permite adaptar el sistema a tipos de análisis que eran imposibles de prever en el instante en que se creó este programa. Se puede pensar en el concepto de extensibilidad total.

5.2 CONCEPTO DE TIPO DE PROBLEMA

Se define tipo de problema como la información que hay que proporcionar al sistema para que este pueda operar con un programa de análisis determinado. O sea, que dado un sistema de preproceso genérico como el desarrollado, y dado un programa de cálculo, se entiende como tipo de problema a toda la información que falta para unir ambos programas.

Figura 67: Esquema de interacción entre preproceso, proceso y tipo de problema.

Esta información consiste en un conjunto de ficheros de texto que describen todos los datos particulares del análisis que el usuario deberá proporcionar posteriormente para cada análisis que vaya a realizar. También habrá una descripción de la manera en que deben proporcionarse esos datos al programa de análisis. Por tanto, la creación de un tipo de problema es un proceso que se realiza de manera completamente externa al sistema de preproceso y que se hace una sola vez. Una vez terminada, el preproceso y el proceso quedan perfectamente unidos y la percepción que recibe el usuario es la de un programa único.

5.3 ABSTRACCIÓN DE DATOS

Cuando pensamos en un análisis por métodos numéricos (elementos finitos, diferencia finitas, etc.), los datos que se necesitan para éste se consideran, en un primer nivel, como entes con significación física. Así, una fuerza o una velocidad inicial o un material tienen un significado físico muy determinado dentro de un análisis.

En un segundo nivel, puede pensarse cada conjunto de estos datos como un conjunto de nombres, números, vectores y matrices que tienen una estructura concreta y distinta para cada agrupación de estos datos. Así, una fuerza puntual se entendería como un vector tridimensional más un punto de aplicación. Un material sería un conjunto de propiedades representables en su mayoría por números, aunque también podría haber propiedades representadas con curvas de variación o con nombres específicos.

El tercer nivel de abstracción consiste en simplificar toda esta información en un conjunto de campos, donde cada campo contiene el nombre del dato a introducir y su valor es siempre representable mediante una cadena de texto. Con este criterio, una fuerza no es más que un conjunto de tres campos, de nombres Componente x, Componente y y Componente z¹ y de valores tres cadenas de texto que en este caso representan números.

En realidad, es necesario introducir unos pocos tipos de campos y mediante éstos será posible representar cualquier información para el análisis. Los tipos de campos son:

1. Campo de valor libre. Se permite al usuario entrar cualquier valor en él. Su uso típico es para entrar valores numéricos o no numéricos para una propiedad.
2. Campos con número fijo de valores. El usuario puede escoger entre un número finito de posibilidades. Su representación gráfica es mediante un menú o unos check boxes².
3. Campo de ejes locales. Permite asociar unos ejes locales a la condición a la que están asociados.
4. Campo de vectores o matrices. Se usa típicamente para representar una propiedad mediante una curva definida por puntos. También es útil para propiedades cuya definición requiere de un número de valores no conocidos a priori.

Con este conjunto de tipos de datos se puede representar virtualmente cualquier tipo de propiedad o información asociada a un análisis.

(¹) O cualquier otros nombres similares.

(²) Se denomina check boxes a un tipo de botones cuyas dos alternativas son estar conectado o desconectado.

5.4 APLICACIÓN A LOS DATOS DEL ANÁLISIS

La separación de los datos del análisis en subgrupos ha sido ya dada en el Apartado 2.5. Cada uno de estos subgrupos consistirá en un conjunto de campos de los anteriormente definidos y en una posible información adicional relativa al grupo en si. A continuación se describe la aplicación de estos criterios a cada uno de los grupos.

5.4.1 Condiciones

Una condición es un conjunto de campos de los anteriormente descritos que contiene una vinculación con las entidades geométricas. Una vez definida, el usuario asignará unos valores particulares de dichos datos a un conjunto de entidades geométricas o de malla. Debe tenerse en cuenta toda la casuística que provoca que estos datos estén finalmente vinculados a un conjunto de nodos, de caras de elementos o de elementos.

NUMBER: 1 CONDITION: Point-Constraints
CONDTYPE: over points
QUESTION: X-Constraint#CB#(0,1)
VALUE: 1
QUESTION: X-Value
VALUE: 0.0
QUESTION: Y-Constraint#CB#(0,1)
VALUE: 1
QUESTION: Y-Value
VALUE: 0.0
QUESTION: Z-Constraint#CB#(0,1)
VALUE: 1
QUESTION: Z-Value
VALUE: 0.0
END CONDITION

Figura 68: Comparación entre los datos que definen una condición en el tipo de problema con la ventana que pide esos datos al usuario.

En la Figura 68 se describen los datos necesarios para definir una condición concreta. En este caso se trata de una prescripción de movimiento nodal con posibilidad de introducir desplazamientos impuestos.

5.4.2 Materiales

Son muy similares al concepto de condición. La diferencia básica consiste en que el tipo de problema proporciona una base de datos de materiales que, posteriormente, el usuario podrá modificar y ampliar.

Los materiales siguen siendo un conjunto de los campos descritos previamente que típicamente representarán a las propiedades de dicho material que tengan interés para un análisis concreto.

5.4.3 Datos generales de problema y de intervalo

Conjunto de campos que para el caso de los datos generales del problema es un conjunto único y para el caso de los datos de intervalo, es un conjunto que se repite para cada intervalo.

Como en todos los casos anteriores, el problemas se reduce a definir cada conjunto como una suma de dichos campos.

5.5 CONEXIÓN PREPROCESO-PROCESO

Hasta ahora, se ha descrito el tratamiento de los datos para el análisis en el interior del sistema de preproceso. Una alternativa consistiría en concluir el problema en este punto y escribir todos estos datos en un formato de fichero neutro. Dada la complejidad y variedad que representa esta información, la solución así escogida exigiría la creación de muy complejos programas de interfase de difícil mantenimiento. Además de esta dificultad, la vertebración de un sistema en múltiples y pequeños componentes conlleva una carga de mantenimiento muy difícil de asumir.

En este trabajo se presenta una forma alternativa de resolver este problema. Se propone la extensión del tipo de problema para que también incluya una descripción del formato del fichero de entrada de datos para el cálculo. Se parte de la hipótesis de que el programa de análisis recibe la información de entrada mediante un conjunto de ficheros en los que se incluye toda la información de malla, condiciones de contorno y otros datos.

La manera de hacerlo ha consistido en la creación de un lenguaje de programación o de macros, especializado en esta tarea y que permite describir la forma exacta en que deben escribirse estos datos en el fichero. Los criterios con los que opera este lenguaje se describirán en el siguiente apartado.

5.6 LENGUAJE BAS DE DESCRIPCIÓN DEL FICHERO DE ENTRADA DEL ANÁLISIS

El lenguaje de programación o de macros para definir los datos para el análisis, consiste en un fichero, al que llamaremos fichero BAS, donde se incluye la información que se requiere en el fichero de comunicación. Insertados en medio del texto genérico, que se reproduce íntegramente, se introducen un conjunto de comandos¹ que extraen información de la base de datos y la substituyen en el fichero de salida.

$DIMENSIONES DEL PROBLEMA
DIMENSIONES:  NPNOD= *npoin , NELEM= *nelem , \
NMATS= *nmats , \
NNODE= *nnode , NDIME= *ndime , \ 
NCARG= *GenData(Número_de_cargas) , \
NGDLN= *GenData(Grados_de_libertad)
$

$DIMENSIONES DEL PROBLEMA
DIMENSIONES:  NPNOD=     30 , NELEM=     36 , \
NMATS=      0 , \
NNODE=      3 , NDIME=      2 , \ 
NCARG= 1 , \
NGDLN= 3
$

GEOMETRIA 
$ CONECTIVIDADES ELEMENTALES
$ ELEM. MATER.    SECUENCIA DE CONECTIVIDADES
*loop elems
*elemsnum *elemsmat *elemsConec
*end elems
$ COORDENADAS DE PUNTOS NODALES
$  NODO  COOR.-X   COOR.-Y   COOR.-Z
*loop nodes
*format "%6i%15.5f%15.5f%15.5f"
*NodesNum *NodesCoord
*end
FIN_GEOMETRIA

GEOMETRIA 
$ CONECTIVIDADES ELEMENTALES
$ ELEM. MATER.    SECUENCIA DE CONECTIVIDADES
1      0      1    22     4
2      0     22    21    23
...
36      0     25    30     5
$ COORDENADAS DE PUNTOS NODALES
$  NODO  COOR.-X   COOR.-Y   COOR.-Z
1        -9.94231        2.65128
2        -0.79538        1.50239
...
30        -7.04511        1.39341
FIN_GEOMETRIA

En la Tabla 6, se incluyen como ejemplo algunos comandos que extraen información de propiedades de la malla así como de datos generales del problema. En la Tabla 7, se puede ver otro ejemplo donde se introducen ya estructuras clásicas en la programación como son los bucles. Además de bucles, el lenguaje incluye condicionales, variables y otros tipos de estructuras necesarias para producir la salida adecuada.

Al ser un lenguaje especializado para el cometido de escribir este tipo de información, su programación es mucho más sencilla que la de escribir un programa independiente, que actúe como interfase, en un lenguaje de programación convencional como puede ser FORTRAN o C++.

El lenguaje es interpretado, en el sentido de que cada vez que se va a escribir un fichero de datos para el análisis, se vuelve a leer el fichero BAS y se ejecutan los comandos contenidos en él. La ventaja evidente de esta forma de trabajar reside en que la velocidad de creación y prueba de nuevos tipos de problemas aumenta considerablemente. Se ha comprobado que la velocidad de creación de un nuevo tipo de problema, para un usuario entrenado y para un programa de análisis de dificultad media, puede estimarse en unas pocas horas.

El inconveniente de los lenguajes interpretados, reside en la menor velocidad de ejecución. Ello implica que el tiempo de escritura del fichero de datos es mayor del que tendría una solución especializada. La comprobación práctica, es que el tiempo final para realizar este cometido no representa un problema real para su uso práctico, incluso para modelos de grandes dimensiones.

(¹) Se distingue un comando porque lleva como signo diferenciador el símbolo estrella ($*$) delante suyo.

5.7 EXTENSIONES MEDIANTE LENGUAJES DE SCRIPT

Al tener incluido un lenguaje de script en su interior, el sistema permite la creación de ventanas especializadas y la implementación de funcionalidades nuevas al creador del tipo de problema. Con ello se consigue una alta personalización del sistema que, en su extremo, puede llevar a dar una nueva apariencia al programa.

5.8 APORTACIONES EN LA ADAPTACIÓN AL ANÁLISIS

La idea de la simplificación de los datos en unos pocos tipos básicos, encaja muy bien en la necesidad de adaptación total del sistema. La solución opuesta de tener todos los datos especializados para todos y cada uno de los tipos de análisis impide esta generalización tan beneficiosa.

Finalmente, la solución de insertar un lenguaje especializado y simplificado para escribir los datos en el formato del análisis, es una aportación significativa de esta tesis que cambia el sentido de la adaptación. Del clásico criterio de que el programa de análisis se adaptaba al preproceso, se ha ido al extremo contrario de que el preproceso se convierte en una herramienta básica, a la que todos los programas de análisis se van conectando.

5.9 CONCLUSIONES A LA ADAPTACIÓN

La capacidad de adaptar el sistema de preproceso a virtualmente cualquier clase de análisis representa una importante aportación para todo el colectivo de investigadores que requieren una herramienta de este tipo. Debe tenerse en cuenta que son pocos los programas de cálculo, que por su generalidad y por su uso ampliamente difundido, pueden tener un sistema de pre y postproceso propio de alta calidad. La progresiva especialización de los programas de cálculo exige el disponer de un sistema general de pre y postproceso que pueda adaptarse a sus necesidades. La metodología desarrollada en la tesis y explicada en este capítulo permite llevar a cabo tal cometido con una cantidad de esfuerzo totalmente aceptable.

6 Ejemplos de aplicación

En este capítulo se describirá mediante ejemplos los retos que han permitido al programa de preproceso ir mejorando su calidad y nivel de utilización.

Cada uno de estos ejemplos, que pueden considerarse descriptivos de toda una categoría de modelos, ha supuesto para la tesis la subida de un escalón en el nivel de calidad en las técnicas desarrolladas y en las posibilidades de tratamiento de modelos industriales por el sistema GiD, resultante de este trabajo.

En la descripción unitaria de cada uno de ellos se va a hacer hincapié en las mejoras que ha requerido el modelo para poder ser tratado con eficiencia y robustez.

6.1 Avión caza F16

Este avión fue uno de los primeros casos de importación de modelos externos. La generación de malla de superficie fue usada para análisis de impacto contra una pared. Posibilitó el desarrollo de las NURBS y de la generación de malla sobre ellas. La importación de la geometría se realizó mediante la utilidad de lectura de ficheros batch y el uso de un conversor externo. El modelo está descrito por unas 200 superficies y la malla generada se compone de 30000 triángulos (ver Figura 69).

6.2 Diente de excavadora

Este primer ejemplo (Figura 70), representa la generación de malla sobre un modelo de sólido de un diente de excavadora. Con el se inició la generación de malla de tetraedros sobre formas geométricas no triviales. Las mallas resultantes se componían de 20-40000 tetraedros y los modelos de 30-40 superficies al principio a 300-400 actualmente. Los modelos de fundición completos para la fabricación de estos dientes, como los de las Figuras 71 y 72, pueden llegar a tener varios miles de superficies. Al principio de la tesis representaban un trabajo manual tedioso de preparación de la geometría y de asignación de tamaños de malla. En la actualidad, la importación de la geometría y la generación de malla de tetraedros sobre uno de estos dientes con el sistema GiD, es una operación directa en la que ha desaparecido el trabajo manual. Los análisis realizados sobre ellos han sido tanto de fundición con un modelo termo-mecánico acoplado, como de sólido con contactos entre piezas.

6.3 Presa bóveda

Por los requerimientos del análisis, convenía obtener una malla de hexaedros de 20 nodos sobre el modelo de la presa de la Figura 73. Para ello, se usaron las funcionalidades del sistema GiD en relación a la malla estructurada. En concreto, fue necesario crear funciones que permitieran la creación automática de volúmenes estructurados y la asignación automática y coherente de tamaños y transiciones entre ellos. En la Figura 74, se observa la malla estructurada resultante del modelo. El cálculo se realizó mediante un análisis de sólido con un modelo de degradación de los materiales [25].

6.4 Catedral de Tarazona

El modelo geométrico tridimensional de la catedral de Tarazona se realizó íntegramente en el interior del sistema GiD. Ello puso a prueba las capacidades de modelación para una geometría no sencilla. La malla final se generó con unos 40000 tetraedros. Un requerimiento de este análisis era la minimización del número de elementos con una distribución óptima en las zonas de mayor interés estructural. El cálculo consistió en un análisis estructural con un modelo de daño [25]. En las Figuras 75 y 76, puede verse el modelo geométrico y la malla obtenida.

6.5 Estampación de chapa

Para análisis de estampación de chapa es necesario obtener la malla de los útiles que intervienen en el proceso [26]. Estos útiles representan formas complejas como pueden ser componentes de automóviles. El tratamiento de estos componentes requirió una gran mejora en la importación de ficheros IGES así como en la generación de malla sobre NURBS. Debe tenerse en cuenta que este tipo de piezas contiene radios de acuerdo y otros pequeños detalles que provocan superficies extremadamente complejas con ratios de relación de tamaños de ${\textstyle 1/1000}$ y mayores. Además, con estos modelos se profundizó en la generación de mallas no conexas entre superficies, con malla estructurada automática, con criterios de error cordal y con altas distorsiones con el objetivo último de describir bien la geometría con el mínimo número de elementos.

Los modelos típicos de esta metodología contienen de varios cientos a varios miles de superficies y las mallas obtenidas se sitúan entre 40000 y 400000 elementos triangulares. En la Figura 77, puede verse un ejemplo de tales modelos.

6.6 Transbordador espacial

Para la generación de malla volumétrica alrededor del transbordador espacial de la Figura 78, se importó la geometría mediante una superficie-malla, se creó la caja contenedora y se remalló el conjunto. Ello supuso la implementación de la generación de malla de superficie sobre una malla previa. El análisis posterior consistió en un cálculo fluido dinámico (CFD) con superficie libre que representaba la forma de las olas al zambullirse el transbordador en el mar [27].

6.7 Gran Telescopio de Canarias

El Gran Telescopio de Canarias es un edificio abierto, en el que interesaba modelizar el movimiento del aire en su interior. Para ello se necesitaba la descripción geométrica del edificio más la descripción topográfica del terreno circundante. Mediante la importación y mejora posterior de la geometría del edificio, más la importación de la topografía del terreno se creó el dominio de análisis. La topografía se incorporó mediante cuadrículas de puntos. Se requirió para ello desarrollar la lectura de ficheros de batch y mejorar la generación automática de superficies.

La malla de tetraedros del dominio, del orden de 700000 tetraedros, requirió un incremento de nivel de calidad de la generación volumétrica para poder generar malla en las ventanas de unión interior-exterior, alrededor del telescopio en si, y realizar transiciones de tamaño hacia el contorno exterior. En las Figuras 79 y 80, puede apreciarse parte de la geometría del modelo.

6.8 Análisis de velas de barco

El análisis de velas de barco se realiza con un modelo acoplado CFD-estructural con interpolación de presiones del fluido hacia el modelo de sólido. La creación del modelo geométrico requería la creación de superficies suaves a partir de un conjunto de líneas paralelas y el duplicado de los nodos alrededor de las velas para el cálculo fluido dinámico. Esta última opción se ha visto facilitada por el criterio de asimilar topología geométrica a topología de malla. De todas maneras, es necesario subdividir el volumen en varios subvolúmenes para evitar la problemática de elementos cruzados sobre las velas, ya que estas últimas son, en el fondo, unos agujeros colapsados en el interior del dominio.

En las Figuras 81 y 82, puede verse el modelo y el contorno de la malla generada.

6.9 Análisis de un barco

El modelo geométrico del casco de barcos se realiza normalmente con un número pequeño de NURBS, donde cada una de ellas tiene una forma muy sofisticada. En este caso, la generación de malla convencional sobre un dominio alterado era incapaz de producir elementos de calidad sobre el casco del barco. Con la generación directa de triángulos en el espacio, se ha podido reproducir bien la forma con una buena relación de aspecto. Las Figuras 83 y 84, muestran uno de los modelos y su malla. Las mallas para análisis hidrodinámico sobre cascos de barco oscilan entre 300000 y un millón de tetraedros.

6.10 Fundición de aluminio

En la pieza de la Figura 85, se puede ver un modelo para análisis de procesos de fundición de aluminio. La característica de más interés a nivel de generación de malla es que es posible mallar correctamente un volumen importado mediante fichero IGES con cientos de superficies y complejidad muy alta. Las mallas para estos análisis rondan los cientos de miles de elementos. Estos análisis tienen la complejidad adicional de que hay que crear contactos pieza-molde, tanto mecánicos como térmicos.

6.11 Importación de la geometría de una fragata

A nivel de comprobación de la calidad de la importación de geometría externa, así como de su tratamiento y visualización en el interior del sistema, se presenta este ejemplo de una fragata con altísimo nivel de detalle (Figura 87), cuya forma queda descrita con más de 10000 superficies.

Figura 69: Malla de un avión caza F16.

Figura 70: Malla de una pieza de diente de excavación, recortada para permitir ver los tetraedros interiores.

Figura 71: Geometría de un conjunto de piezas de diente de excavación. Ésta es la forma en que se construyen, por fundición.

Figura 72: Visualización foto-realística de la geometría anterior.

Figura 73: Geometría de una presa de bóveda.

Figura 74: Malla de la presa anterior generada como malla estructurada de elementos hexaédricos.

Figura 75: Geometría de la catedral de Tarazona.

Figura 76: Malla de la catedral.

Figura 77: Análisis de procesos de estampación de chapa.

Figura 78: Transbordador espacial entrando en el mar.

Figura 79: Malla del edificio del telescopio de Gran Canarias para realizar un análisis fluido dinámico (CFD).

Figura 80: Malla y resultados del campo de velocidades del aire alrededor del telescopio con parte del terreno circundante.

Figura 81: Geometría de una vela de barco deportivo más el dominio de análisis.

Figura 82: Malla de la vela anterior.

Figura 83: Geometría de un barco deportivo más el dominio de análisis.

Figura 84: Malla del barco anterior.

Figura 85: Malla de una pieza creada para la fundición de aluminio.

Figura 86: Detalle de la malla anterior.

Figura 87: Representación geométrica de un barco importado desde IGES.

7 Conclusiones y líneas futuras de investigación

7.1 PRINCIPALES APORTACIONES DE LA TESIS

Cada uno de los capítulos concluye con un apartado destinado a especificar y distinguir lo que se considera aportación de la tesis en contraposición a la algoritmia extraída del estado del arte. En el presente apartado se incluye un extracto de dichas aportaciones.

En el Capítulo 2 se han introducido un conjunto de criterios de diseño del sistema y se han organizado los datos que se requieren para suministrarlos a los códigos de análisis, con el propósito de acercarse al óptimo de facilidad y eficiencia de uso (ver Apartado 2.7).

En el Capítulo 3 se presenta el esquema de estructura jerárquica y las implicaciones que conlleva para generar una malla con los condicionantes que requiere el modelo de análisis. Asimismo, se describen un conjunto de algoritmos destinados a mejorar la calidad de la geometría en cuanto a las necesidades del mallado posterior. También se incluyen otros algoritmos desarrollados para poder actuar sobre la geometría, tanto para crear modelos como para prepararlos para la generación de malla (ver Apartado 3.11).

En el Capítulo 4 se describen un conjunto de algoritmos de generación de malla donde se introducen las aportaciones mediante una necesaria descripción de los métodos generales de mallado. En el Apartado 4.19 se distingue de manera detallada lo que se considera aportación en el conjunto de toda la metodología.

el Capítulo 5 expone los criterios generales y las técnicas que permiten la adaptación del sistema a un código particular de análisis. En concreto, se considera aportación al hecho de suministrar una subdivisión de los datos en unos subconjuntos precisos así como de permitir la configuración desde el exterior mediante un conjunto de macros y de un lenguaje de programación especializado (ver Apartado 5.8).

7.2 CONCLUSIONES GENERALES

En cada uno de los capítulos se ha incluido un apartado de conclusiones que hacía referencia a toda la temática descrita en dicho capítulo. Estas últimas conclusiones generales serán un resumen de aportaciones de todo lo ya descrito, más un conjunto de observaciones que afectan al conjunto de la tesis como unidad.

En primer lugar, para cumplir los objetivos propuestos, era necesario crear capacidades de tratamiento de modelos geométricos complejos. Estas capacidades incluyen tanto la importación a partir de programas externos como una cierta capacidad de modificación de los modelos importados. Además, también debe ser posible la creación dentro del sistema de modelos geométricos hasta una cierta complejidad.

Se ha conseguido este propósito, ya que es habitual con el sistema desarrollado el tratamiento de modelos con varios cientos e incluso miles de superficies NURBS con una relativa comodidad. Se ha dotado al sistema con un conjunto de herramientas que permiten modificar hasta cierto punto esos modelos con el objetivo de convertirlos en aptos para su análisis. Esas mismas herramientas y otras permiten crear modelos de dificultad baja y media sin depender de un programa externo. Es posible crear modelos bi y tridimensionales del orden de dificultad del análisis resistente de una catedral, del análisis aerodinámico de un modelo de un telescopio más el terreno circundante o de las velas de una embarcación recreativa.

En segundo lugar, es necesario corregir la geometría para poder generar una malla de calidad suficiente sobre ella.

La mayoría de los modelos que se importan para crear una malla, han sido diseñados con otros propósitos y por tanto, no optimizados para el análisis. En la tesis, se han creado un conjunto de herramientas que corrigen gran cantidad de estos defectos sin intervención del usuario.

En tercer lugar, se requería la generación de malla de manera eficiente y que se adaptase a gran variedad de necesidades. El resultado ha sido:

• Es posible generar malla estructurada y no estructurada con diversas tipologías de elementos. las velocidades de generación de malla no structurada de triángulos están en el orden de magnitud de 150000 elementos por minuto en un ordenador pentium medio. Las de tetraedros en unos 30000 elementos por minuto para mallas de tamaños típicos.
• Es posible generar malla de superficie sobre piezas industriales, como componentes de automóviles, que mantengan criterios de calidad geométricos y que minimicen el número de elementos finales. El número de elementos generados típico para esas mallas oscila entre 50000-200000 elementos.
• La generación de malla de tetraedros se realiza en problemas de fluidodinámica para análisis de barcos, velas, telescopios y en modelos termomecánicos para estudios de piezas industriales de fundición. Los tamaños típicos de malla oscilan entre 200000-1000000 elementos.
• Existe la posibilidad de creación de malla especializada para contactos, juntas, simetrías rotacionales, multi-volúmenes superpuestos, etc.

Finalemente, se requería que el programa fuese capaz de tratar todos los datos de diferentes tipologías de programas de análisis. Al mismo tiempo, el coste de adaptación entre programas debía ser bajo y realizable de manera exterior al sistema de preproceso.

El resultado ha sido que actualmente existe una gran cantidad de códigos tanto académicos como comerciales que se han adaptado completamente al sistema GiD desarrollado en la tesis. Después de la adaptación, la percepción que recibe el usuario es la de un programa único. Ejemplos de estas adaptaciones son:

• Programas de análisis de sólidos tanto estáticos como dinámicos y con no linealidad tanto geométrica como de materiales. Calsef [28], COMET [25].
• Programas de estampación industrial y de impacto. Stampack [26], Simpack [29].
• Programas termo-mecánicos de fundición. COMET [25].
• Programas térmicos. CALTEP [30].
• Análisis magnetoestáticos. EMANT [31].
• Programas de CFD. FLASH [32], PHANTOM [27].
• Programas de análsis de rádar mediante técnicas de ray tracing modificado. GRECO.
• Análisis acoplado fluido-térmico. PHANTOM [27].
• Programas de diseño de velas por métodos numéricos. SAILS.
• Programas de mecánica del suelo. BRITE-EURAM.
• Otros.

7.3 LÍNEAS FUTURAS DE INVESTIGACIÓN

A continuación se incluye una lista de posibles temas de interés para investigación futura dentro del ámbito de trabajo planteado:

1. Desarrollar nuevas técnicas de mejora automática de la geometría con el objetivo de generar una malla sobre ella. Estos algoritmos deben ser, en su mayor parte, de ejecución sin la intervención del usuario.
2. Algunas técnicas de modelación geométrica, como la intersección de superficies o la creación de radios de acuerdo espaciales, por su propia complejidad, pueden considerarse como un tema aún no completamente resuelto. Esta línea de trabajo profundizaría en el tratamiento de entidades con grados elevados de distorsión.
3. La generación de malla no estructurada de hexaedros es un tema no resuelto que sería altamente necesario para algunos tipos de modelos que no pueden operar correctamente con tetraedros.
4. La generación de capas límite para análisis mediante Navier-Stokes provoca una serie de problemáticas que deben ser resueltas de una manera totalmente general. Entre éstas se destaca la necesidad de generar una malla estructurada próxima a las superficies que definen el modelo de análisis para capturar la capa límite, que esté conectada a una malla no estructurada que discretice el dominio exterior.
5. El uso de superficies-malla debe extenderse para mallas no conformes y para conjuntos de puntos en el espacio. En concreto, a partir de una nube de puntos puede interpolarse un conjunto de superficies NURBS o una malla de triángulos que puede usarse más tarde para mallar sobre ella. Esta capacidad facilitaría el modelado de objetos descritos por puntos a partir de maquinaria de reconocimiento de formas tridimensionales.
6. La creación de geometría a partir de tomografías médicas facilitaría la entrada de la ingeniería en el mundo de la medicina.
7. Las técnicas de adaptación del sistema a un programa de cálculo pueden mejorarse para permitir mayor versatilidad de uso por parte del usuario. Puede pensarse en conceptos tales como:
• Permitir la ayuda interactiva para las propiedades a asignar.
• Criterios lógicos de interacción de campos. Unos campos pueden forzar valores determinados a otros campos.
• La claridad se consigue con la separación de un problema complejo en varios de más sencillos. Deben darse opciones de subdivisión de los campos en diversos componentes.

8 Descripción matemática de las NURBS

8.1 INTRODUCCIÓN

En los siguientes apartados se va a describir de manera general la formulación matemática de las líneas y superficies NURBS¹. Se pretende dar una idea general de los principales algoritmos y formulaciones necesarios para modelar con estas entidades y generar malla sobre ellas, sin acompañarlos con extensas demostraciones teóricas o explicaciones detalladas. Para una descripción más profunda del tema ver [1]. Muchas de las ideas de este apéndice se han obtenido a partir de la referencia anterior.

Previamente, se definirán algunas propiedades básicas del espacio euclídeo que servirán de introducción a la nomenclatura y propiedades de las curvas.

(¹) Acrónimo de Non Uniform Rational B-Splines.

8.2 ESPACIO EUCLÍDEO

Dado un conjunto de puntos y un conjunto de vectores, diremos que los puntos pertenecen al espacio euclídeo tridimensional ${\textstyle E^{3}}$ y que los vectores pertenecen a ${\textstyle R^{3}}$.

Dados los puntos ${\textstyle {\vec {P_{0}}}}$ y ${\textstyle {\vec {P_{1}}}}$, la diferencia entre ambos da como resultado el vector ${\textstyle {\vec {v}}}$. La substracción de dos elementos de ${\textstyle E^{3}}$ da como resultado un elemento de ${\textstyle R^{3}}$.

${\vec {v}}={\vec {P_{0}}}-{\vec {P_{1}}}\qquad {\vec {P_{0}}},{\vec {P_{1}}}\in E^{3},\quad {\vec {v}}\in R^{3}$

En general la operación de suma para puntos no está bien definida. Sólo es válida la combinación baricéntrica, que consiste en la suma de un conjunto de puntos afectados cada uno de ellos por un coeficiente o peso. La suma de los pesos debe ser igual a la unidad.

${\vec {P}}=\sum _{i=0}^{n}\alpha _{i}{\vec {P_{i}}}\qquad \left|{\begin{array}{l}{\vec {P_{i}}}\in E^{3}\\\sum _{i=0}^{n}\alpha _{i}=1\end{array}}\right.$

Esta operación es válida debido a que puede considerarse como la suma de un punto con un vector:

${\vec {P}}={\vec {P_{0}}}+\sum _{i=1}^{n}\alpha _{i}({\vec {P_{i}}}-{\vec {P_{0}}})$

Una combinación baricéntrica es convexa si todos los coeficientes ${\textstyle \alpha _{i}}$ son mayores o iguales a cero. El lugar geométrico de todos los puntos resultantes de una combinación de este tipo se llama polígono convexo. Una propiedad de este polígono es que la recta que une dos puntos cualesquiera de su interior, estará también totalmente contenida en su interior. En la Figura 88, se aprecia una recta cualquiera que une dos puntos extremos.

Figura 88: Polígono convexo del conjunto de puntos.

En general, todas las curvas y superficies que se describirán en los apartados siguientes, se construyen a partir de combinaciones baricéntricas de puntos. De aquí la importancia de esta introducción.

8.3 DIFERENTES TIPOS DE LÍNEAS

Para describir el tipo de líneas denominadas NURBS en toda su generalidad, será interesante introducir previamente un conjunto de curvas más simples que nos servirán como introducción a algunas de sus propiedades. En concreto, el orden natural de introducción de las curvas será:

1. Curvas de Bézier. Tipo de curvas que se definen mediante un polígono de control con un número fijo de puntos en función del grado.
2. Splines. Extensión de las anteriores a partir del concepto de curva definida a trozos. Incorporan los knots¹.
3. B-splines. Otro nombre dado a las Splines cuando se representan mediante el polígono de Boor.
4. Cónicas. Curvas con término racional que describen al conjunto de: parábolas, hipérbolas, elipses, círculos, etc.
5. Líneas NURBS. B-splines con término racional. O sea, que combinan las propiedades de las B-splines con las de las cónicas.

Todas las curvas que trataremos serán paramétricas tridimensionales, siendo el espacio paramétrico en todos los casos el ${\textstyle [0,1]}$. La introducción de esta restricción en el espacio del parámetro no resta generalidad a la formulación, pues se podrían obtener ecuaciones equivalentes en el dominio genérico ${\textstyle [a,b]}$, con la simple aplicación de una transformación del parámetro ${\textstyle t}$ tal que ${\textstyle t'=a+(b-a)t}$.

(¹) Se llama knots al conjunto de parámetros $K_{i}$ pertenecientes a $R$ tal que $K_{i}\in [0,1]$, que forman una serie de valores posibles del parámetro de la curva y que sirven para su definición. Una posible traducción sería nodos. En todo caso, en este documento mantendremos el vocablo original por razones de uso común.

8.4 CURVAS DE BÉZIER

Las curvas de Bézier son una forma alternativa de representación de una curva polinómica. En lugar de la representación tradicional de un polinomio como: ${\textstyle {\vec {P}}={\vec {P}}_{0}+{\vec {P}}_{1}t+\dots +{\vec {P}}_{n}t^{n}}$, se escoge otra representación en la que los puntos ${\textstyle {\vec {P}}^{i}}$ adquieren sentido geométrico. A continuación se introduce, de manera constructiva, la definición de una curva de Bézier.

Figura 89: Curvas de Bézier de diversos grados.

Dados los puntos ${\textstyle {\vec {P_{0}}}}$ y ${\textstyle {\vec {P_{1}}}}$, todos los puntos ${\textstyle {\vec {P}}}$ pertenecientes a ${\textstyle E^{3}}$ de la forma

${\vec {P}}={\vec {P}}(t)=(1-t){\vec {P_{0}}}+t{\vec {P_{1}}}\qquad t\in R$

se denominan la línea recta a través de ${\textstyle {\vec {P_{0}}}}$ y ${\textstyle {\vec {P_{1}}}}$. Para ${\textstyle 0\leq t\leq 1}$, el punto ${\textstyle {\vec {P}}}$ estará entre los dos puntos. Por tanto podemos considerar a ${\textstyle {\vec {P}}(t)}$ como una función paramétrica de ${\textstyle t\in [0,1]}$ hacia un segmento de línea recta en ${\textstyle E^{3}}$. Nótese que esta función se obtiene como combinación baricéntrica de dos puntos.

Si realizamos la misma construcción para un grado más, o sea, creamos una nueva función de ${\textstyle t}$ tal que sea combinación baricéntrica de dos combinaciones baricéntricas lineales, se obtiene:

${\vec {P_{0}^{1}}}(t)$ $=(1-t){\vec {P_{0}}}+t{\vec {P_{1}}}$ ${\vec {P_{1}^{1}}}(t)$ $=(1-t){\vec {P_{1}}}+t{\vec {P_{2}}}$ Insertamos ambas funciones en una nueva combinación baricéntrica, ${\vec {P_{0}^{2}}}(t)$ $=(1-t){\vec {P_{0}^{1}}}+t{\vec {P_{1}^{1}}}$ y desarrollamos, obteniendo la parábola: ${\vec {P_{0}^{2}}}(t)$ $=(1-t)^{2}{\vec {P_{0}}}+2t(1-t){\vec {P_{1}}}+t^{2}{\vec {P_{2}}}$

Por tanto, obtenemos una nueva función paramétrica, llamada parábola, que es una aplicación del intervalo ${\textstyle [0,1]}$ a una curva de grado dos en el espacio. Algunas de sus propiedades pueden deducirse directamente de las propiedades de la combinación baricéntrica y son:

1. La curva pasa por el punto inicial y por el final.
2. La tangente de la curva en el punto inicial y en el final pasan por el punto central.
3. La curva es continua y derivable, de grado ${\textstyle C^{1}}$ en todo el dominio.
4. La curva es plana y está incluida en el polígono convexo que forman sus puntos.

Si generamos la construcción anterior para grado ${\textstyle n}$, obtendremos:

${\vec {P}}_{i}^{r}(t)=(1-t){\vec {P}}_{i}^{r-1}(t)+t{\vec {P}}_{i+1}^{r-1}(t)\qquad {\hbox{con: }}\quad \left|{\begin{array}{rl}r=&1,\dots ,n\\i=&0,\dots ,n-r\end{array}}\right.$

(113)

Definimos la curva de Bézier de grado ${\textstyle n}$, ${\textstyle {\vec {P}}^{n}}$ como ${\textstyle {\vec {P}}_{0}^{n}(t)}$ en (113). Vemos que esta curva queda perfectamente definida mediante su grado ${\textstyle n}$, y su polígono de control, que son los puntos ${\textstyle {\vec {P}}_{0},\dots ,{\vec {P}}_{n}}$. Esta curva cumple las siguientes propiedades:

1. La curva pasa por el punto inicial y por el final.
2. La tangente de la curva en el ${\textstyle {\vec {P}}_{0}}$ y ${\textstyle {\vec {P}}_{n}}$ pasan, respectivamente, por los puntos ${\textstyle {\vec {P}}_{1}}$ y ${\textstyle {\vec {P}}_{n-1}}$ si no coinciden con los anteriores y, en general, por los más cercanos a ellos que no sean coincidentes.
3. La curva es continua y derivable, de grado ${\textstyle C^{n-1}}$ en todo el dominio.
4. La curva está incluida en el polígono convexo que forma su polígono de control.
5. De las propiedades anteriores se deduce que la curva se acerca suavemente a los puntos interiores del polígono de control.
6. Se deduce también que si el polígono de control reside en un plano, la curva estará contenida en el mismo plano.

En la Figura 89, pueden observarse las propiedades de interpolación, continuidad y tangencias para curvas de Bézier de diversos grados.

El método usado para introducir y definir a las curvas de Bézier es al mismo tiempo, una función recursiva que sirve para evaluarlas. Existen otras formas de definir y evaluar este tipo de curvas que se describirán en apartados posteriores.

8.4.1 Curvas de Bézier expresadas en términos de los polinomios de Bernstein

Se definen los polinomios de Bernstein ${\textstyle B_{i}^{n}(t)}$ como:

$B_{i}^{n}(t)$ $={\binom {n}{i}}t^{i}(1-t)^{n-i}$ Donde se recuerda que: ${\binom {n}{i}}$ $=\left\{{\begin{array}{ll}{\frac {n!}{i!(n-i)!}}&{\hbox{si }}0\leq i\leq n\\0&{\hbox{en caso contrario. }}\end{array}}\right.$

Los polinomios de Bernstein cumplen las siguientes propiedades:

• Recursividad.

$B_{i}^{n}(t)=(1-t)B_{i}^{n-1}(t)+tB_{i-1}^{n-1}(t)\quad {\hbox{con: }}\left|{\begin{array}{l}B_{0}^{0}(t)\equiv 1\\B_{j}^{n}(t)\equiv 0\quad {\hbox{para }}j\not \in [0,\dots ,n]\end{array}}\right.$

• El conjunto forma una partición de la unidad.

$\sum _{i=0}^{n}B_{i}^{n}(t)\equiv 1$

Dadas estas definiciones, podemos expresar una curva de Bézier ${\textstyle {\vec {P}}^{n}(t)}$ como:

${\vec {P}}^{n}(t)={\vec {P}}_{0}^{n}(t)=\sum _{i=0}^{n}{\vec {P}}_{i}B_{i}^{n}(t).$

8.4.2 Diversas propiedades de las curvas de Bézier

8.4.2.1 Derivadas

Si definimos el operador diferencias en avance ${\textstyle \Delta }$ como ${\textstyle \Delta {\vec {P}}_{i}={\vec {P}}_{i+1}-{\vec {P}}_{i}}$, obtenemos que la derivada de una curva de Bézier es:

${\frac {d}{dt}}{\vec {P}}^{n}(t)=n\sum _{i=0}^{n}\left[B_{i-1}^{n-1}(t)-B_{i}^{n-1}(t)\right]{\vec {P}}_{i}=n\sum _{i=0}^{n-1}\Delta {\vec {P}}_{i}B_{i}^{n-1}(t)\qquad {\hbox{con }}\Delta {\vec {P}}_{i}\in R^{3}$

(115)

Para las derivadas de orden superior definimos el operador diferencias en avance iterativo ${\textstyle \Delta ^{r}}$ como:

$\Delta ^{r}{\vec {P}}_{i}=\Delta ^{r-1}{\vec {P}}_{i+1}-\Delta ^{r-1}{\vec {P}}_{i}=\sum _{j=0}^{r}{\binom {r}{j}}(-1)^{r-j}{\vec {P}}_{i+j}$

(116)

A partir de (115) y (116) se deduce que la derivada de grado ${\textstyle r}$ de ${\textstyle {\vec {P}}^{n}(t)}$ es:

${\frac {d^{r}}{dt^{r}}}{\vec {P}}^{n}(t)={\frac {n!}{(n-r)!}}\sum _{i=0}^{n-r}\Delta ^{r}{\vec {P}}_{i}B_{i}^{n-r}(t)$

(117)

Operando, es posible deducir una expresión recursiva para las derivadas que mantendría una relación directa con (113). Para obtener sus expresiones ver [1].

8.4.2.2 Subdivisión

Dada una curva ${\textstyle {\vec {P}}(t)}$ con ${\textstyle t\in [0,1]}$, queremos encontrar una curva ${\textstyle {\vec {P}}'(t)}$ tal que ${\textstyle {\vec {P}}(t)\equiv {\vec {P}}'(s)}$ para ${\textstyle s\in [0,s_{1}]}$ siendo ${\textstyle s_{1}}$ un valor del parámetro situado en el interior del intervalo ${\textstyle [0,1]}$. O sea, queremos subdividir la curva ${\textstyle {\vec {P}}(t)}$ y obtener una nueva curva que sea equivalente a ella en un subconjunto del intervalo del parámetro. Si ${\textstyle {\vec {P}}_{0},\dots ,{\vec {P}}_{n}}$ es el polígono de control de ${\textstyle {\vec {P}}(t)}$ y ${\textstyle {\vec {P}}'_{0},\dots ,{\vec {P}}'_{n}}$ es el polígono de control de ${\textstyle {\vec {P}}'(s)}$, entonces debe cumplirse que:

${\vec {P}}'_{i}={\vec {P}}_{0}^{i}(s_{1})$

(118)

Donde el superíndice ${\textstyle i}$ se refiere a que son los puntos del polinomio de grado ${\textstyle i}$ calculado a partir de ${\textstyle {\vec {P}}_{0},\dots ,{\vec {P}}_{n}}$ que es el polinomio de grado 0.

En la Figura 90, se observa el nuevo polígono de control calculado a partir del inicial y del parámetro ${\textstyle s_{1}}$.

Figura 90: Subdivisión de curvas de Bézier. A partir de la curva de los ${\vec {P}}_{i}$ definida en $[0,1]$, obtenemos la curva definida por los ${\vec {P}}'_{i}$, definida en $[0,s_{1}]$.

8.4.2.3 La forma matricial

Otra forma de expresar las curvas de Bézier es mediante matrices en la siguiente forma:

${\vec {P}}^{n}(t)=\sum _{i=0}^{n}{\vec {P}}_{i}B_{i}^{n}(t)=\left[{\vec {P}}_{0},\dots ,{\vec {P}}_{n}]\right]\left[{\begin{array}{c}B_{0}^{n}(t)\\\vdots \\B_{n}^{n}(t)\end{array}}\right]$ $=\left[{\vec {P}}_{0},\dots ,{\vec {P}}_{n}]\right]\left[{\begin{array}{ccc}m_{00}&\cdots &m_{0n}\\\vdots &&\vdots \\m_{n0}&\cdots &m_{nn}\end{array}}\right]\left[{\begin{array}{c}t^{0}\\\vdots \\t^{n}\end{array}}\right]$ donde los coeficientes $m_{ij}$ valen: $m_{ij}$ $=(-1)^{j-i}{\binom {n}{j}}{\binom {j}{i}}$

8.4.2.4 Esquema de Horner

Una forma más eficiente de evaluar una curva de Bézier es mediante el esquema de Horner, que se basa en anidar las operaciones unas dentro de otras. Para el caso de una curva de grado cúbico tendremos:

${\vec {P}}^{3}(t)=\left\{\left[{\binom {3}{0}}s{\vec {P}}_{0}+{\binom {3}{1}}t{\vec {P}}_{1}\right]s+{\binom {3}{2}}t^{2}{\vec {P}}_{2}\right\}s+{\binom {3}{3}}t^{3}{\vec {P}}_{3}$

con ${\textstyle s=1-t}$ y teniendo en cuenta que ${\textstyle {\binom {n}{i}}={\frac {n-i+1}{i}}{\binom {n}{i-1}}}$ para ${\textstyle i>0}$.

8.4.2.5 Incremento del grado de la curva

Dada una curva de Bézier ${\textstyle {\vec {P}}^{n}(t)}$ de grado ${\textstyle n}$, queremos obtener otra curva ${\textstyle {\vec {P}}^{'n+1}(t)}$ de grado ${\textstyle n+1}$ tal que ${\textstyle {\vec {P}}^{n}(t)\equiv {\vec {P}}^{'n+1}(t)}$. Mediante interpolación lineal a trozos y operando se obtiene el nuevo polígono de control ${\textstyle {\vec {P}}'_{0},\dots ,{\vec {P}}'_{n+1}}$ que debe cumplir:

${\vec {P}}'_{i}={\frac {i}{n+1}}{\vec {P}}_{i-1}+(1-{\frac {i}{n+1}}){\vec {P}}_{i}\qquad {\hbox{para }}i=0,\dots ,n+1$

(125)

Si generalizamos la ecuación anterior para un incremento de grado ${\textstyle r}$, se deduce fácilmente que:

${\vec {P}}_{i}^{'(r)}=\sum _{j=0}^{n}{\vec {P}}_{j}{\binom {n}{j}}{\frac {\binom {r}{i-j}}{\binom {n+r}{i}}}$

(126)

La demostración detallada de este proceso puede encontrarse en [1].

8.4.2.6 Integrales de los polinomios de Bernstein y de las curvas de Bézier

Se puede comprobar que:

$\int _{0}^{t}B_{i}^{n}(x)\,dx$ $={\frac {1}{n+1}}\sum _{j=i+1}^{n+1}B_{j}^{n+1}(t)$ $\int _{0}^{1}B_{i}^{n}(x)\,dx$ $={\frac {1}{n+1}}$ $\int _{0}^{t}{\vec {P}}^{n}(x)\,dx$ $={\frac {t}{n+1}}\sum _{i=0}^{n}{\vec {P}}_{0}^{i}(t)$

8.5 SPLINES

Tal como se ha descrito en la sección 8.4, el grado de una curva de Bézier determina el número de puntos que contiene el polígono de control. Debido a esto, cuando se necesita aumentar el número de puntos para tener más flexibilidad en la definición de la curva, se debe aumentar el grado. A efectos prácticos y de coste computacional, grados excesivamente altos (${\textstyle n>10}$) no son aconsejables. La extensión natural de estas curvas para obtener otras más flexibles, es la de las curvas polinómicas a trozos. Se trata de definir una nueva curva que sea la unión de varias curvas polinómicas simples (como las curvas de Bézier) y que cumplan una serie de requisitos de continuidad en la unión entre cada dos de estas curvas polinómicas. Esta nueva curva se llamará spline.

Figura 91: Una curva spline se forma mediante la unión de sucesivas curvas de Bézier cumpliendo unos requisitos de continuidad en los puntos de unión.

Si tenemos ${\textstyle L}$ curvas de Bézier unidas de manera consecutiva y cada una de ellas definida en el intervalo ${\textstyle [0,1]}$, tal como se muestra en la Figura 91, podemos realizar una transformación local del parámetro tal que queden definidas en ${\textstyle [a_{i},b_{i}]}$. Así obtenemos una curva global ${\textstyle {\vec {P}}_{G}(t)}$ como unión de todas ellas y cuyo parámetro estará definido en ${\textstyle [a,b]}$. Se tiene que cumplir que:

$a=a_{0};\qquad b=b_{L-1};\qquad a_{i}=b_{i-1};\qquad {\hbox{y }}\qquad b-a=\sum _{i=0}^{L-1}b_{i}-a_{i}.$

Si llamamos ${\textstyle u_{0},\dots ,u_{L}}$ a ${\textstyle a_{0},a_{1},\dots ,a_{L-1},b_{L-1}}$, La curva ${\textstyle {\vec {P}}_{G}(t)}$ quedará perfectamente definida con los ${\textstyle u_{i}}$ y con el conjunto de los polígonos de control de las curvas de Bézier:

${\vec {P}}_{00},\dots ,{\vec {P}}_{0n},\dots ,{\vec {P}}_{L-1,0},\dots ,{\vec {P}}_{L-1,n}\qquad {\hbox{para }}L{\hbox{ curvas de grado }}n{\hbox{. }}$

Por continuidad ${\textstyle C^{0}}$ de las curvas se comprueba que ${\textstyle {\vec {P}}_{in}={\vec {P}}_{i+1,0}}$. El conjunto de los ${\textstyle u_{i}}$ se llama la secuencia de knots. En la nota 8.3, se comenta sobre la posible traducción de este término.

Nota: El tamaño de cada intervalo ${\textstyle [a_{i},b_{i}]}$ en relación al tamaño total, es una elección arbitraria que cambiará la parametrización de la curva. Los efectos de esta elección se comentarán más adelante.

En general, además de la continuidad ${\textstyle C^{0}}$ de las curvas, se deseará que tengan continuidad en las derivadas. Más concretamente, nos interesará que el punto de unión entre dos segmentos de curva tenga el mismo grado de continuidad que las curvas mismas.

Continuidad $\mathbf {C^{1}}$

Para que la curva sea continua ${\textstyle C^{1}}$ en todo el dominio de definición, imponemos igualdad de derivadas en la unión de dos segmentos. Si llamamos ${\textstyle \Delta _{i}}$ a ${\textstyle u_{i+1}-u_{i}}$, ${\textstyle t}$ al parámetro global de la curva y ${\textstyle t'_{i}}$ a los parámetros locales de cada segmento, se deduce que:

${\frac {d}{dt}}{\vec {P}}_{G}(t)={\frac {1}{\Delta _{0}}}{\frac {d}{dt'_{0}}}{\vec {P}}_{0}(t'_{0})={\frac {1}{\Delta _{1}}}{\frac {d}{dt'_{1}}}{\vec {P}}_{1}(t'_{1})\quad {\hbox{como: }}\left|{\begin{array}{l}{\frac {d}{dt'_{0}}}{\vec {P}}_{0}(1)=n\Delta {\vec {P}}_{n-1}\\{\frac {d}{dt'_{1}}}{\vec {P}}_{1}(0)=n\Delta {\vec {P}}_{n}\end{array}}\right.\quad {\hbox{entonces: }}$

$\Delta _{1}\Delta {\vec {P}}_{n-1}=\Delta _{0}\Delta {\vec {P}}_{n}$

(127)

Se observa del resultado que para que la curva sea derivable ${\textstyle C^{1}}$ en todo el dominio, es condición necesaria pero no suficiente que los tres puntos del polígono de control alrededor de un punto de unión deben ser colineales. Las relaciones de distancia a las que deben estar entre ellos depende de la parametrización. Por tanto, dada una curva formada por segmentos de Bézier y con los puntos colineales 3 a 3 alrededor de los puntos de unión, será ${\textstyle C^{1}}$ en todo el dominio o no, en función de la parametrización escogida.

Continuidad $\mathbf {C^{2}}$

Suponiendo la curva ${\textstyle {\vec {P}}_{G}(t)}$ continua ${\textstyle C^{1}}$, la condición de continuidad ${\textstyle C^{2}}$ en el punto ${\textstyle {\vec {P}}_{n}}$ implica que debe existir un punto ${\textstyle {\vec {P}}_{S}}$ que cumpla:

${\vec {P}}_{n-1}$ $=(1-\alpha ){\vec {P}}_{n-2}+\alpha {\vec {P}}_{S}$ (128) ${\vec {P}}_{n+1}$ $=(1-\alpha ){\vec {P}}_{S}+\alpha {\vec {P}}_{n+2}$ (129)

Con ${\textstyle \alpha =\Delta _{0}/(t_{2}-t_{0})}$. En la Figura 92, puede verse el criterio de continuidad ${\textstyle C^{2}}$. La deducción de estas fórmulas se obtiene fácilmente a partir de la imposición de igualdad de derivasdas segundas.

Figura 92: La condición de continuidad $C^{2}$ en el punto ${\vec {P}}_{n}$ implica la unicidad del punto ${\vec {P}}_{S}$

A partir de la aplicación de estos conceptos, se puede comprobar que la forma geométrica de la curva queda perfectamente definida con los polígonos de control de los segmentos de Bézier. Por otro lado, la continuidad de la curva variará en función de la parametrización escogida. Por tanto, con una misma forma geométrica se pueden obtener curvas con diferentes grados de continuidad en función del tipo de parametrización escogida. Debe notarse que una transformación del tipo ${\textstyle t'=\alpha {+\beta }t}$ para toda la secuencia de knots, no cambia ni la forma de la curva ni su grado de continuidad.

8.6 B-SPLINES

Las curvas B-spline son curvas spline representadas y definidas mediante el polígono B-spline y la secuencia de knots. Seguidamente se va a deducir la forma del polígono B-spline para curvas de grado ${\textstyle C^{1}}$ y ${\textstyle C^{2}}$ de manera constructiva a partir de las curvas splines y, posteriormente, se va a dar la definición genérica de B-spline para cualquier grado.

8.6.1 Construcción de la B-spline cuadrática ($\mathbf {C^{1}}$)8.6.1 Construcción de la B-spline cuadrática ( C 1 {\displaystyle \mathbf {C^{1}} } )

Supongamos que tenemos una curva spline cuadrática, continua ${\textstyle C^{1}}$ en todo su dominio y compuesta por ${\textstyle L}$ segmentos de Bézier. La curva queda definida mediante los polígonos de control ${\textstyle {\vec {P}}_{00},{\vec {P}}_{01},{\vec {P}}_{02},\dots ,{\vec {P}}_{L0},{\vec {P}}_{L1},{\vec {P}}_{L2}}$ y la secuencia de knots. Para que se cumpla el criterio de continuidad ${\textstyle C^{1}}$ en los puntos ${\textstyle {\vec {P}}_{2i}}$, se requiere que:

${\vec {P}}_{2i}={\frac {\Delta _{i}}{\Delta _{i-1}-\Delta _{i}}}{\vec {P}}_{2i-1}+{\frac {\Delta _{i-1}}{\Delta _{i-1}-\Delta _{i}}}{\vec {P}}_{2i+1}$

(130)

Por tanto, el mínimo número de puntos para que quede perfectamente definida esta curva es el polígono ${\textstyle {\vec {P}}_{00},{\vec {P}}_{01},{\vec {P}}_{03},\dots ,{\vec {P}}_{L-1,1},{\vec {P}}_{L1},{\vec {P}}_{L2}}$. A esta secuencia de puntos se le llama el polígono B-spline o polígono de Boor. Una curva B-spline cuadrática quedará definida por el polígono B-spline y por la secuencia de knots. En la Figura 93, puede apreciarse dicha relación.

Figura 93: Relación entre los puntos de control de la B-spline y de las curvas de Bézier que la forman, para el caso cuadrático.

8.6.2 Construcción de la B-spline cúbica ($\mathbf {C^{2}}$)8.6.2 Construcción de la B-spline cúbica ( C 2 {\displaystyle \mathbf {C^{2}} } )

Para que una curva spline cúbica sea continua ${\textstyle C^{2}}$ en todos los puntos ${\textstyle {\vec {P}}_{3i},\,i\geq 0}$, debe cumplirse el criterio de continuidad ${\textstyle C^{1}}$ y además, debe existir un conjunto de puntos ${\textstyle {\vec {P}}_{Sj},\,j\geq 1}$ que cumplan:

${\vec {P}}_{3i}={\frac {\Delta _{i}}{\Delta _{i-1}-\Delta _{i}}}{\vec {P}}_{3i-1}+{\frac {\Delta _{i-1}}{\Delta _{i-1}-\Delta _{i}}}{\vec {P}}_{3i+1}$ (131) ${\vec {P}}_{3i-2}={\frac {\Delta _{i-1}+\Delta _{i}}{\Delta }}{\vec {P}}_{S,i+1}+{\frac {\Delta _{i-2}}{\Delta }}{\vec {P}}_{S,i+2}$ (132) ${\vec {P}}_{3i-1}={\frac {\Delta _{i}}{\Delta }}{\vec {P}}_{S,i+1}+{\frac {\Delta _{i-2}+\Delta _{i-1}}{\Delta }}{\vec {P}}_{S,i+2}$ (133)

con ${\textstyle \Delta =\Delta _{i-2}+\Delta _{i-1}+\Delta _{i}}$. Esta imposición viene dada de que se fuerza igualdad de derivadas segundas en los puntos. Cerca de los extremos de la línea, usaremos ${\textstyle {\vec {P}}_{S,1}={\vec {P}}_{0}}$ y ${\textstyle {\vec {P}}_{S,2}={\vec {P}}_{1}}$. Entonces se obtiene que:

${\vec {P}}_{2}={\frac {\Delta _{1}}{\Delta _{0}+\Delta _{1}}}{\vec {P}}_{S,2}+{\frac {\Delta _{0}}{\Delta _{0}+\Delta _{1}}}{\vec {P}}_{S,3}$

Para los puntos finales de la línea se procede de igual manera. El conjunto de los puntos ${\textstyle {\vec {P}}_{S,1},\dots ,{\vec {P}}_{S,L+3}}$ es el polígono B-spline de la curva. Se observa que el número de puntos de control de la B-spline cúbica ${\textstyle {\mathit {NP}}_{s}={\mathit {NP}}_{b}/3+3}$ siendo ${\textstyle NP_{b}}$ el número de puntos del polígono spline.

Se ha obtenido, tanto para el caso cuadrático como para el cúbico, una formulación que nos permite calcular los polígonos de control de los segmentos Bézier a partir del polígono B-spline. De esta manera podemos evaluar la curva para cualquier ${\textstyle t\in [0,1]}$. En el siguiente apartado se generaliza el concepto de B-spline para cualquier grado y se describen técnicas de evaluación de las curvas directamente, sin tener que acceder a los segmentos de Bézier.

Se puede observar en la Figura 94, la relación entre los polígonos de control de la curva de Bézier con el polígono B-spline.

Figura 94: Relación del polígono de control de una B-spline cúbica con las curvas de Bézier asociadas.

8.6.3 Definición de las curvas B-spline

Una curva B-spline viene definida por su grado, su polígono de control y su secuencia de knots.

Como una B-spline es una curva formada por diversos segmentos polinómicos, el grado de la línea será el mismo grado que tendrá cada uno de los segmentos.

Definición: El orden de una curva es el grado más uno.

El polígono de control es un conjunto de puntos en el espacio donde algunos de ellos pueden estar repetidos. La curva pasará por el primer y por el último punto y se acercará, de manera suave, a todos los otros puntos.

La secuencia de knots, es un conjunto de números reales en forma de lista no decreciente. Normalmente se definen en el intervalo ${\textstyle [0,1]}$ y, por tanto, el primero valdrá 0 y el último 1. Si un mismo valor está repetido ${\textstyle r}$ veces, se dice que tiene una multiplicidad ${\textstyle r}$. Diferentes multiplicidades implicarán diferentes propiedades en relación al punto de control que corresponde al knot. En concreto, una multiplicidad igual al orden supondrá que la curva interpole al punto correspondiente. El primer knot de la lista y el último tendrán multiplicidad igual al orden para asegurar la interpolación de los puntos extremos.

El número de puntos y el de knots deben cumplir que:

$N_{k}=N_{p}+o\qquad {\hbox{con: }}\left|{\begin{array}{rl}N_{k}=&{\hbox{Número de ''knots'' }}\\N_{p}=&{\hbox{Número de puntos }}\\o=&{\hbox{orden }}\end{array}}\right.$

Para evaluar una curva de este tipo se pueden usar varios algoritmos. Uno de ellos proviene de la inserción de knots y se enuncia de la siguiente manera: Dada una curva definida por su grado ${\textstyle n}$, polígono de control ${\textstyle {\vec {P}}_{1},\dots ,{\vec {P}}_{L}}$ y knots ${\textstyle u_{0},\dots ,u_{L+n}}$, para evaluarla en ${\textstyle u}$ tendremos que encontrar ${\textstyle u_{I}}$ tal que: ${\textstyle u\in [u_{I},u_{I+1}]}$. El valor de la curva en ${\textstyle u}$ será:

${\vec {P}}_{i}^{k}(u)={\frac {u_{i+n-k}-u}{u_{i+n-k}-u_{i-1}}}{\vec {P}}_{i-1}^{k-1}(u)+{\frac {u-u_{i-1}}{u_{i+n-k}-u_{i-1}}}{\vec {P}}_{i}^{k-1}(u)\quad {\hbox{para: }}\left|{\begin{array}{l}k=1,\dots ,n-r\\i=I-n+k+1,\dots ,I+1\\{\vec {P}}_{i}^{0}(u)={\vec {P}}_{i}\end{array}}\right.$ entonces: ${\vec {s}}(u)={\vec {P}}_{I+1}^{n-r}(u)$

Donde ${\textstyle r}$ es la multiplicidad del knot ${\textstyle u_{I}=u}$. Si no existe ese knot, entonces ${\textstyle r=0}$. Se deduce de esta ecuación que la evaluación de la curva pasa por calcular recursivamente una serie de interpolaciones lineales.

Nota: Obsérvese que en la notación usada, la lista de puntos empieza con el subíndice 1 y la lista de knots empieza con el 0. Si se quisiera que ambas listas empezasen por el mismo subíndice, deberían modificarse ligeramente las fórmulas.

8.6.4 Ejemplo de evaluación de una línea B-spline cuadrática

Para ilustrar los puntos que intervienen en la evaluación de la curva para un valor del parámetro dado, se va a dar un pequeño ejemplo de una curva B-spline cuadrática ${\textstyle (n=2)}$, definida mediante un polígono de control de 4 puntos ¹. El polígono de control y la secuencia de knots se pueden ver en la Figura 95.

Figura 95: Ejemplo de B-spline cuadrático con 4 puntos de control y evaluada en $u$.

Suponiendo que ${\textstyle u\in [u_{3},u_{4})}$, entonces ${\textstyle I=3}$ y la curva ${\textstyle {\vec {s}}(u)={\vec {P}}_{4}^{2}}$. El cálculo será el siguiente:

${\vec {P}}_{3}^{1}=(1-\alpha ){\vec {P}}_{2}+\alpha {\vec {P}}_{3}\qquad {\hbox{con: }}\alpha ={\frac {u-u_{2}}{u_{4}-u_{2}}}$ (136) ${\vec {P}}_{4}^{1}=(1-\alpha ){\vec {P}}_{3}+\alpha {\vec {P}}_{4}\qquad {\hbox{con: }}\alpha ={\frac {u-u_{3}}{u_{5}-u_{3}}}$ (137) ${\vec {P}}_{4}^{2}=(1-\alpha ){\vec {P}}_{3}^{1}+\alpha {\vec {P}}_{4}^{1}\qquad {\hbox{con: }}\alpha ={\frac {u-u_{3}}{u_{4}-u_{3}}}$ (138)

Se comprueba que para una curva cuadrática, la función recursiva se evalúa en dos niveles y que los puntos de control que afectan al resultado son:

${\begin{array}{c}{\vec {P}}_{2}\\{\vec {P}}_{3}\\{\vec {P}}_{4}\end{array}}{\begin{array}{c}{\vec {P}}_{3}^{1}\\{\vec {P}}_{4}^{1}\end{array}}{\vec {P}}_{4}^{2}$

y que los knots afectados son ${\textstyle u_{2},u_{3},u_{4},u_{5}}$. A medida que aumenta el grado de la curva, aumentan los niveles de recursividad y también el número de puntos de control que intervienen en la solución.

(¹) El mínimo número de puntos que puede tener el polígono de control de una B-spline cuadrática es de 3. En este caso, sería equivalente a una curva de Bézier.

8.6.5 Parametrizaciones

En diferentes casos nos encontraremos con que no disponemos de la secuencia de knots y que hay que calcularla. Un caso típico sería cuando un usuario crea interactivamente una línea dando los puntos de control o modificándolos. El programa debe calcular los knots por si solo. Otro caso en el que debe calcularse esta lista es en la interpolación de puntos mediante B-splines.

El problema se enuncia como: Dado el polígono de control de la curva ${\textstyle {\vec {P}}_{1},\dots ,{\vec {P}}_{L}}$ y su grado ${\textstyle n}$, calcular una posible secuencia de knots, ${\textstyle u_{0},\dots ,u_{L+n}}$ tal que la curva resultante sea la deseada¹. Este cálculo de parámetros puede realizarse con diversas técnicas, de entre las que se comentan dos:

Parametrización constante Todos los incrementos del parámetro son iguales. Es la más sencilla. Si la separación entre puntos es muy distinta de unos a otros, puede dar avances muy poco constantes.

Ejemplo: ${\textstyle 0,0,0,0,0.2,0.4,0.6,0.8,1,1,1,1}$ para ${\textstyle n=3}$ y ${\textstyle L=8}$.

Parametrización chord length El incremento relativo del parámetro es igual a la relación de la distancia entre puntos de control². Así, ${\textstyle u_{0}=0}$,${\textstyle \Delta u_{i}=\left\|\Delta {\vec {P}}_{i}\right\|/u_{L+n}}$. Esta parametrización introduce información geométrica dentro de la lista y da, en general, mucho mejores resultados. En particular, el avance en la curva al incrementar el parámetro es mucho más constante.

Ejemplo: para ${\textstyle n=3}$ y ${\textstyle L=6}$ se tendría ${\textstyle 0,0,0,0,d_{1}/d_{3},d_{2}/d_{3},1,1,1,1}$ con ${\textstyle d_{1}=\left\|{\vec {P}}_{3}-{\vec {P}}_{1}\right\|}$, ${\textstyle d_{2}=d_{1}+\left\|{\vec {P}}_{4}-{\vec {P}}_{3}\right\|}$ y ${\textstyle d_{3}=d_{2}+\left\|{\vec {P}}_{6}-{\vec {P}}_{4}\right\|}$.

Para interpolación de puntos mediante B-splines (ver 8.6.10), los algoritmos se aplicarán a los puntos a interpolar en lugar de a los puntos de control. Esto es debido a que se necesita conocer la lista de knots, para poder calcular los puntos de control.

(¹) Normalmente, consideraremos una curva mejor, si los cambios de curvatura son suaves. También se tendrá en cuenta que el parámetro de la curva se aproxime al parámetro arco. O sea, que el avance en la curva sea lo más constante posible.

(²) El segundo y penúltimo puntos no se tienen en cuenta para este cálculo de distancias.

8.6.6 La base B-spline

De manera similar a como definíamos los polinomios de Bernstein en la sección 8.4.1, también se puede definir una base de polinomios para el espacio de las B-spline. Esta base nos permite encontrar otra forma de evaluar las curvas. Si consideramos a las ${\textstyle N_{I}^{n}(u)}$ como funciones recursivas definidas de la siguiente manera:

$N_{i}^{0}(u)=\left\{{\begin{array}{l}1\qquad {\hbox{si }}u_{i-1}\leq u<u_{i}\\0\qquad {\hbox{en caso contrario. }}\end{array}}\right.$ $N_{i}^{n}(u)={\frac {u-u_{i-1}}{u_{i-n-1}-u_{i-1}}}N_{i}^{n-1}(u)+{\frac {u_{i+n}-u}{u_{i+n}-u_{i}}}N_{i+1}^{n-1}(u)$ entonces: ${\vec {s}}(u)=\sum _{i=1}^{n}{\vec {P}}_{i}N_{i}^{n}$

(140)

Entre las dos técnicas de evaluación de la curva, esta última requiere menos operaciones que la anterior debido a que la recursividad afecta solamente a términos escalares. Por otro lado, la técnica primera es más estable numéricamente.

8.6.7 Derivada de una curva B-spline

La derivada de una curva B-spline puede deducirse a partir de la expresión de la curva y su valor es:

${\frac {d}{du}}{\vec {s}}(u)$ $={\frac {n}{u_{I+1}-u_{I}}}\Delta {\vec {P}}_{I}^{n-1}$ ${\frac {d}{du}}{\vec {s}}(u)$ $=n\sum _{i=1}^{L+n-1}{\frac {\Delta {\vec {P}}_{i-1}}{u_{n+i-1}-u_{i-1}}}N_{i}^{n-1}(u)$

Para calcular la derivada expresada en función de los polinomios de la base B-spline, se usará la expresión:

${\frac {d}{du}}N_{i}^{n}(u)={\frac {n}{u_{n+i-1}-u_{i-1}}}N_{i}^{n-1}(u)-{\frac {n}{u_{n+i}-u_{i}}}N_{i+1}^{n-1}(u)$

8.6.8 Efectos de la multiplicidad de knots y puntos sobre la continuidad de la curva

En la Tabla 8, se describen los efectos en la continuidad que provoca la repetición de knots o puntos de control para una B-spline de grado 3. Debe tenerse en cuenta que La notación ${\textstyle G^{2}}$, se refiere a que la curva tiene continuidad geométrica de grado 2 pero no necesariamente continuidad en la derivada segunda. En [2], se describen gráficamente los efectos sobre la curva de las diferentes multiplicidades.

Tabla. 8 Efectos de la multiplicidad para una B-spline cúbica.

Multiplicidad	Puntos de control múltiples	knots múltiples
1	$C^{2}$,${\textstyle G^{2}}$	$C^{2}$,${\textstyle G^{2}}$
2	$C^{2}$,${\textstyle G^{1}}$	$C^{1}$,${\textstyle G^{1}}$
3	$C^{2}$,${\textstyle G^{0}}$ La curva interpola el punto de control triple. Los segmentos en cada lado del punto son líneas rectas.	$C^{0}$,${\textstyle G^{0}}$ La curva interpola el punto de control. Se puede controlar la forma de los segmentos a cada lado del punto.
4	$C^{2}$,${\textstyle G^{0}}$ La curva interpola el punto de control cuádruple. Los segmentos en cada lado del punto son líneas rectas y interpolan a los otros puntos de control.	Hay discontinuidad en la curva. La curva acaba en un punto de control y continúa en el siguiente. Se puede controlar la forma de los segmentos a cada lado del punto.

8.6.9 Representación gráfica de una B-spline

Cuando se desea dibujar en una pantalla gráfica o en un fichero un conjunto de entidades geométricas, la primitiva de dibujo de que se dispone en la librería gráfica es la de pintado de líneas rectas. Por tanto, para representar una curva ésta debe descomponerse en un conjunto elevado de segmentos rectos. La obtención de estos segmentos puede realizarse de dos maneras:

• Evaluar la curva a pequeños intervalos del parámetro ${\textstyle t}$.
• Insertar un conjunto de knots mediante la ecuación (136), y dibujar el polígono de control resultante. El nuevo polígono de control representará a la curva con mayor precisión a medida que se vayan insertando más knots. El número de knots insertados en cada zona puede depender de la curvatura de la curva en esa misma zona.

El segundo método permite dibujar la curva con menos operaciones para igual precisión.

8.6.10 Interpolación mediante B-splines

Existen muchas técnicas para definir una curva que interpole a un conjunto de puntos dados. En este caso se va a realizar una interpolación mediante B-splines cúbicas. La formulación del problema viene dada por: dados los puntos ${\textstyle {\vec {x}}_{1},\dots ,{\vec {x}}_{L}}$, encontrar una curva B-spline de grado 3 y ${\textstyle L+2}$ puntos de control que interpole a los puntos dados. En la Figura 96, puede observarse la relación entre ambos conjuntos de puntos.

Figura 96: Obtención de los puntos de control a partir de un conjunto de puntos interpolados con una B-spline cúbica.

Para ello, debemos crear primero la lista de knots. Las técnicas de creación serán las descritas en la sección 8.6.5.

A partir de las relaciones entre los sucesivos polígonos de Bézier con el polígono de control B-spline (ver [[#8.6.1 Construcción de la B-spline cuadrática ($\mathbf {C^{1}}$)|8.6.1]] y [[#8.6.2 Construcción de la B-spline cúbica ($\mathbf {C^{2}}$)|8.6.2]]), se pueden deducir las siguientes relaciones:

$\Delta _{i}$ $=\Delta u_{i}$ $\alpha _{i}$ $={\frac {\Delta _{i}^{2}}{\Delta _{i-2}+\Delta _{i-1}+\Delta _{i}}}$ $\beta _{i}$ $={\frac {\Delta _{i}(\Delta _{i-2}+\Delta _{i-1})}{\Delta _{i-2}+\Delta _{i-1}+\Delta _{i}}}+{\frac {\Delta _{i-1}(\Delta _{i}+\Delta _{i+1})}{\Delta _{i-1}+\Delta _{i}+\Delta _{i+1}}}$ $\gamma _{i}$ $={\frac {\Delta _{i-1}^{2}}{\Delta _{i-1}+\Delta _{i}+\Delta _{i+1}}}$

$\left[{\begin{array}{ccccccc}1\\\alpha _{1}&\beta _{1}&\gamma _{1}\\&&&\ddots \\&&&&\alpha _{L-1}&\beta _{L-1}&\gamma _{L-1}\\&&&&&&1\end{array}}\right]\left[{\begin{array}{c}{\vec {P}}_{2}\\{\vec {P}}_{3}\\\vdots \\{\vec {P}}_{L}\\{\vec {P}}_{L+1}\end{array}}\right]=\left[{\begin{array}{c}{\vec {r}}_{a}\\(\Delta _{0}+\Delta _{1}){\vec {x}}_{2}\\\vdots \\(\Delta _{L-2}+\Delta _{L-1}){\vec {x}}_{L-1}\\{\vec {r}}_{b}\end{array}}\right]$

${\vec {P}}_{1}={\vec {x}}_{1}\qquad {\vec {P}}_{L+2}={\vec {x}}_{L}$

Donde ${\textstyle {\vec {r}}_{a}}$ y ${\textstyle {\vec {r}}_{b}}$, que son el segundo y penúltimo punto, deben darse como una información adicional. Una posibilidad es que el usuario los introduzca como una indicación de las tangentes a la curva en los puntos extremos, otra posibilidad es autocalcularlos. En este caso existen diversas técnicas, de la que se comenta una de ellas que es la de los puntos finales de Bessel.

La técnica de los puntos finales de Bessel consiste en calcular la parábola que pasa por los 3 puntos extremos de un lado de la curva y considerar que el punto a calcular debe estar en el interior de esta parábola. Normalmente se escoge un punto convencional de la parábola como por ejemplo ${\textstyle t=0.4}$. Se repetiría dos veces el proceso, una para cada lado.

Otra técnica de cálculo de los puntos incógnitas es la de imponer curvatura nula en los extremos. A este criterio de cálculo se le denomina Condición natural de puntos finales.

El cálculo de la curva interpolante se reduce pues a calcular la lista de knots, calcular los dos puntos de definición de tangentes (si no están ya dados), y resolver un sistema de ecuaciones tridiagonal. Más información sobre interpolación geométrica en [33].

La interpolación de un conjunto de puntos con una curva cerrada implica un ligero cambio en el sistema, ya que para asegurar continuidad ${\textstyle C^{2}}$ entre los puntos inicial y final, hay que incluir algunos términos más en la matriz, dejando esta de ser tridiagonal.

8.7 CURVAS RACIONALES. LAS CÓNICAS

Figura 97: Proyección de una curva polinómica cuadrática en el espacio $E^{3}$ al plano $z=\omega =1$. La curva resultante es una cónica en $E^{2}$.

Existen varias maneras de definir una cónica. Una de las posibles definiciones es:

Definición Una cónica en el espacio ${\textstyle E^{2}}$ es la proyección de una parábola en ${\textstyle E^{3}}$ sobre un plano.

Si proyectamos en el plano ${\textstyle z=1}$ con centro de proyección el origen de coordenadas, a partir de un punto ${\textstyle \left[{\begin{array}{ccc}x&y&z\end{array}}\right]^{T}}$, se obtiene la proyección ${\textstyle \left[{\begin{array}{ccc}x/z&y/z&1\end{array}}\right]^{T}}$. De esta manera se comprueba que una cónica, expresada como una curva de Bézier cuadrática, tiene la forma:

${\vec {s}}(t)={\frac {\omega _{0}{\vec {P}}_{0}B_{0}^{2}(t)+\omega _{1}{\vec {P}}_{1}B_{1}^{2}(t)+\omega _{2}{\vec {P}}_{2}B_{2}^{2}(t)}{\omega _{0}B_{0}^{2}(t)+\omega _{1}B_{1}^{2}(t)+\omega _{2}B_{2}^{2}(t)}}$

donde los ${\textstyle \omega _{i}}$ se llaman los pesos de los puntos de la curva.

Otra forma de expresar el mismo problema es considerar que el punto ${\textstyle \left[{\begin{array}{ccc}x&y&z\end{array}}\right]^{T}}$, tiene como coordenadas homogéneas al punto ${\textstyle \left[{\begin{array}{cccc}\omega x&\omega y&\omega z&\omega \end{array}}\right]^{T}}$, que representa a todos los puntos del hiperespacio ${\textstyle E^{4}}$ que proyectan en él. En este caso:

${\vec {s}}(t)=\sum _{i=0}^{2}{\vec {P}}_{Hi}B_{i}^{2}(t)\qquad {\hbox{con }}{\vec {P}}_{Hi}\in E^{4}$

La relación de un punto ${\textstyle {\vec {P}}}$ con su punto homogéneo asociado ${\textstyle \left[{\begin{array}{cc}\omega {\vec {P}}&\omega \end{array}}\right]^{T}}$, puede observarse en la Figura 97.

8.8 LÍNEAS NURBS

Una línea NURBS es una B-spline con término racional. Las expresiones para evaluar una NURBS serán equivalentes a las de evaluación de las B-splines, añadiendo los pesos ${\textstyle \omega _{i}}$ con su evaluación, también recursiva. En el caso de operar con los puntos en coordenadas homogéneas (en ${\textstyle E^{4}}$), las fórmulas son idénticas. En esta sección se darán las expresiones a aplicar usando puntos ${\textstyle {\vec {P}}\in E^{3}}$.

8.8.1 Evaluación recursiva

Esta fórmula es equivalente a la de la ecuación 136, añadiendo las componentes del término racional. El cálculo de los pesos ${\textstyle \omega _{i}^{k}}$ se realiza también de manera recursiva. Para evaluar la curva en ${\textstyle u}$ se debe encontrar ${\textstyle u_{I}}$ tal que ${\textstyle u\in [u_{I},u_{I+1}]}$,

${\vec {P}}_{i}^{k}(u)={\frac {(1-\alpha _{i}^{k})\omega _{i-1}^{k-1}{\vec {P}}_{i-1}^{k-1}(u)+\alpha _{i}^{k}\omega _{i}^{k-1}{\vec {P}}_{i}^{k-1}(u)}{\omega _{i}^{k}}}\quad {\hbox{para: }}\left|{\begin{array}{l}k=1,\dots ,n-r\\i=I-n+k+1,\dots ,I+1\\{\vec {P}}_{i}^{0}(u)={\vec {P}}_{i}\\\omega _{i}^{0}=\omega _{i}\end{array}}\right.$ (144) $\alpha _{i}^{k}={\frac {u-u_{i-1}}{u_{i+n-k}-u_{i-1}}}$ $\omega _{i}^{k}=(1-\alpha _{i}^{k})\omega _{i-1}^{k-1}+\alpha _{i}^{k}\omega _{i}^{k-1}$ entonces: ${\vec {s}}(u)={\vec {P}}_{I+1}^{n-r}(u)$

8.8.2 Evaluación mediante la base B-spline

De manera completamente equivalente a la ecuación 140, sin más que añadir el término racional obtenemos la expresión de la línea NURBS en función de la base B-spline.

${\vec {s}}(u)={\frac {\sum _{i=1}^{n}\omega _{i}{\vec {P}}_{i}N_{i}^{n}}{\sum _{i=1}^{n}\omega _{i}N_{i}^{n}}}$

8.8.3 Derivada de una línea NURBS

La derivada de una línea NURBS puede expresarse como:

${\frac {d}{du}}{\vec {s}}(u)={\frac {\beta (\omega _{I+1}^{n-1}{\vec {P}}_{I+1}^{n-1}-\omega _{I}^{n-1}{\vec {P}}_{I}^{n-1})-{\frac {\beta (\omega _{I+1}^{n-1}-\omega _{I}^{n-1})}{\omega _{I+1}^{n}}}{\vec {P}}_{I+1}^{n}}{\omega _{I+1}^{n}}}\qquad {\hbox{con: }}\beta ={\frac {n}{u_{I+1}-u_{I}}}$

8.8.4 Definición de la circunferencia como línea NURBS

Para ilustrar la definición de una curva racional mediante líneas NURBS, se ofrece el ejemplo de la Figura 98. Los valores de las coordenadas están dados para el círculo unitario centrado en el origen. Cualquier otro círculo se obtiene mediante traslaciones y escalados. Una elipse tendría una definición equivalente, sin más que transformar el cuadrado que forman los puntos de control en un rectángulo.

boxed figure

Figura 98: En este ejemplo se detallan los puntos de control y los knots necesarios para definir un círculo unitario centrado en el origen.

plain figure

8.8.5 Influencia de los pesos en la curva

El cambio de los pesos de los puntos de control de una NURBS afectarán, en general, a la forma de la curva. En concreto, el incremento del peso de un punto tenderá a atraer la curva hacia ese punto. En la Figura fig:influencia-pesos-NURBS, puede verse la diferencia en el cambio de la forma de la curva según se mueva el punto de control o se incremente el peso de ese mismo punto. Esta propiedad de los pesos puede usarse para modelar la forma de la curva.

Si todos los pesos tienen el mismo valor, la forma de la curva es equivalente a la curva no racional. Esto limita el interés de los pesos para modelar porque, si se varían varios pesos de varios puntos, la curva tiende a recuperar su forma original. Se parte de la restricción de que todos los pesos deben ser positivos para evitar singularidades.

Se dice que la curva está en la forma estándar si los pesos de los dos puntos extremos son iguales a uno. Nos interesará, en general, que la curva esté en esta forma para poder realizar algunas operaciones como creación de superficies a partir de su contorno. Toda curva se puede convertir a su forma estándar mediante:

${\hat {\omega }}_{i}=\left({\sqrt[{n}]{\frac {\omega _{1}}{\omega _{n+1}}}}\right)^{i-1}\quad i=1,\dots ,n+1$ (145) ${\hat {\omega }}_{i}={\frac {1}{\omega _{1}}}\omega _{i}\quad i=1,\dots ,n+1$ (146)

La ecuación (145) está descrita para curvas de Bézier racionales pero puede extenderse a NURBS. La ecuación (146) aplica la propiedad de que un escalado uniforme de los pesos no varía la forma de la curva. Estas transformaciones de los pesos influirán en la parametrización.

8.9 SUPERFICIES NURBS

Si realizamos la extensión de las expresiones de las líneas NURBS en dos direcciones, obtenemos las superficies NURBS. Éstas son el producto tensorial de dos líneas.

Figura 99: En esta figura se puede observar una superficie NURBS con su polígono de control.

Las superficies quedarán perfectamente definidas mediante mediante el grado y la lista de knots en cada una de las direcciones ${\textstyle u}$ y ${\textstyle v}$. Tanto el grado como el número de knots, no tienen que coincidir para ${\textstyle u}$ y ${\textstyle v}$. Los puntos de control serán una cuadrícula en el espacio que en general tendrá diferente número de puntos según cada una de las direcciones.

La expresión general para evaluar una superficie NURBS será:

${\vec {s}}(u,v)={\frac {\sum _{i}\sum _{j}\omega _{ij}{\vec {P}}_{ij}N_{i}^{m}(u)N_{j}^{n}(u)}{\sum _{i}\sum _{j}\omega _{ij}N_{i}^{m}(u)N_{j}^{n}(u)}}$

La mayoría de las operaciones a realizar sobre estas superficies, incluida la evaluación, pueden realizarse primero sobre un sentido, obtener como resultado una línea y aplicar la operación sobre el otro sentido. De esta manera reducimos el problema de operar sobre una superficie al problema de operar sobre un conjunto de líneas.

Ejemplo de evaluación: Dada una superficie con ${\textstyle m}$ puntos en dirección ${\textstyle u}$ y ${\textstyle n}$ puntos en dirección ${\textstyle v}$, que deseamos evaluarla en ${\textstyle (u_{a},v_{a})}$, tenemos que evaluar las ${\textstyle m}$ curvas formadas con los puntos de la cuadrícula de control tomados como columnas y evaluarlas en ${\textstyle v=v_{a}}$. El resultado será un conjunto de ${\textstyle m}$ puntos que formarán una curva que hay que evaluar para ${\textstyle u=u_{a}}$.

8.9.1 Cálculo de la normal

Figura 100: Cuando el lado de una superficie NURBS degenera a un punto, el cálculo de la normal no puede realizarse mediante las derivadas.

El cálculo de la normal se realizará a partir del producto vectorial de las derivadas. Para obtener la normal a la superficie en ${\textstyle (u,v)}$ se calculará de la siguiente manera:

${\vec {N}}={\frac {{\frac {\delta {\vec {S}}(u,v)}{\delta u}}\wedge {\frac {\delta {\vec {S}}(u,v)}{\delta v}}}{\left\|{\frac {\delta {\vec {S}}(u,v)}{\delta u}}\wedge {\frac {\delta {\vec {S}}(u,v)}{\delta v}}\right\|}}$

En el caso de superficies con alguno de sus lados colapsados a un punto, la expresión anterior daría normal 0 al evaluarla en el contorno colapsado. Por tanto, no sería válida. Para estas situaciones, se puede usar la siguiente fórmula para evaluar la normal en el contorno:

${\vec {N}}(0,0)={\frac {\Delta _{u}{\vec {P}}_{11}\wedge \Delta _{v}{\vec {P}}_{11}}{\left\|\Delta _{u}{\vec {P}}_{11}\wedge \Delta _{v}{\vec {P}}_{11}\right\|}}$

(147)

En la ecuación (147), debe hacerse un control, caso por caso, de los puntos que quedan colapsados y modificar ligeramente los puntos a tener en cuenta para el cálculo. En la Figura 100, puede observarse un caso de superficie degenerada.

8.9.2 Interpolación a partir del contorno

Dadas un conjunto cerrado de líneas en el espacio, queremos encontrar una superficie, cuyo contorno sean las líneas dadas, y que interpole de manera suave a esas líneas. El proceso consiste en reducir el conjunto de líneas a uno equivalente formado por cuatro líneas NURBS. Estas las asociamos dos a dos mediante el criterio de líneas opuestas. Cada uno de estos grupos debe tener el mismo grado y la misma lista de knots. Ello se consigue mediante técnicas de elevación de grado y mediante inserción de knots (ecuación 144).

Al acabar este proceso, de la cuadrícula de puntos de control de la superficie conocemos todos los puntos del contorno y conocemos los grados y las listas de knots. Únicamente hay que calcular los puntos interiores de la cuadrícula. Para ello se pueden usar varios algoritmos interpolantes. Se presenta a continuación algunos de ellos.

8.9.2.1 Superficie Coon bilineal

Dado el contorno de una superficie, representado por ${\textstyle {\vec {x}}(u,v)}$ con ${\textstyle u,v\in [0,1]}$, la superficie Coon bilineal se crea mediante la composición de las dos superficies regladas definidas tomando las curvas opuestas dos a dos, y un término de orden superior para compatibilizarlas. Su definición es la siguiente:

${\vec {x}}(u,v)=\left[{\begin{array}{cc}1-u&u\end{array}}\right]\left[{\begin{array}{c}{\vec {x}}(0,v)\\{\vec {x}}(1,v)\end{array}}\right]+\left[{\begin{array}{cc}{\vec {x}}(u,0)&{\vec {x}}(u,1)\end{array}}\right]\left[{\begin{array}{c}1-v\\v\end{array}}\right]$ $-\left[{\begin{array}{cc}1-u&u\end{array}}\right]\left[{\begin{array}{cc}{\vec {x}}(0,0)&{\vec {x}}(0,1)\\{\vec {x}}(1,0)&{\vec {x}}(1,1)\end{array}}\right]\left[{\begin{array}{c}1-v\\v\end{array}}\right]$

Si consideramos el polígono de control de una superficie NURBS como una superficie Coon, y definimos unos valores ${\textstyle u_{i},v_{i}}$ para cada punto de control, podemos calcular los puntos de control interiores mediante esta formulación.

El cálculo de los ${\textstyle u_{i},v_{i}}$ debe realizarse con algún criterio parecido al de chord length (ver sección 8.6.5) porque un espaciamiento uniforme de los parámetros puede producir malos resultados en superficies distorsionadas.

Nota: Los ${\textstyle u_{i},v_{i}}$ referenciados, se usan únicamente para el cálculo de los puntos interiores del polígono de control y, por tanto, no tienen relación alguna con la lista de knots.

8.9.2.2 Superficie Coon parcialmente bicúbica

Las superficies bilineales tienen el problema de que las tangentes en el contorno se calculan con información que esta alejada de él. Como consecuencia, dos superficies cuyas líneas son continuas ${\textstyle C^{1}}$, tienen discontinuidad ${\textstyle C^{1}}$ y ${\textstyle G^{1}}$. En la Figura 101, se observa la discontinuidad entre ambas superficies.

Figura 101: Las dos superficies se han creado a partir del contorno de líneas. Puede observarse una pequeña discontinuidad $G^{1}$ entre ambas.

Una superficie que se comporta mejor ante este efecto es la parcialmente bicúbica, que es una modificación de la anterior que introduce los polinomios cúbicos de Hermite. El uso de estos polinomios procede de la interpolación mediante una curva de Bézier cúbica, de dos puntos y las dos derivadas en ellos. Su expresión es:

$H_{0}^{3}(t)$ $=B_{0}^{3}(t)+B_{1}^{3}(t)$ $H_{0}^{3}(t)$ $={\frac {1}{3}}B_{1}^{3}(t)$ $H_{0}^{3}(t)$ $=-{\frac {1}{3}}B_{2}^{3}(t)$ $H_{0}^{3}(t)$ $=B_{2}^{3}(t)+B_{3}^{3}(t)$

Donde los ${\textstyle B_{i}^{j}}$ son los polinomios de Bernstein definidos en 8.4.1. La expresión de esta superficie será:

${\vec {x}}(u,v)=\left[{\begin{array}{cc}H_{0}^{3}&H_{3}^{3}\end{array}}\right]\left[{\begin{array}{c}{\vec {x}}(0,v)\\{\vec {x}}(1,v)\end{array}}\right]+\left[{\begin{array}{cc}{\vec {x}}(u,0)&{\vec {x}}(u,1)\end{array}}\right]\left[{\begin{array}{c}H_{0}^{3}\\H_{3}^{3}\end{array}}\right]$ $-\left[{\begin{array}{cc}H_{0}^{3}&H_{3}^{3}\end{array}}\right]\left[{\begin{array}{cc}{\vec {x}}(0,0)&{\vec {x}}(0,1)\\{\vec {x}}(1,0)&{\vec {x}}(1,1)\end{array}}\right]\left[{\begin{array}{c}H_{0}^{3}\\H_{3}^{3}\end{array}}\right]$

8.10 SUPERFICIES NURBS TRIMADAS

Las superficies NURBS trimadas¹ se definen mediante una superficie NURBS como base y una o varias curvas de corte. Según el sentido de la curva, ésta será un contorno exterior o un agujero. Típicamente, se obtienen a partir de la intersección de dos superficies NURBS normales. Las líneas de intersección definirán a las líneas de corte.

Figura 102: La figura de la izquierda es una superficie NURBS. La de la derecha es otra superficie NURBS con la misma superficie de base y con unas líneas de recorte formadas por la intersección con el prisma dibujado.

En las Figuras 102 y 103, se observa que a partir de una intersección de una superficie NURBS con otras superficies se puede obtener tanto la superficie externa a las líneas de corte como la interna. En este último caso, nos quedaría la superficie original con un agujero.

Figura 103: En este caso, el recorte provocado por la intersección crea un agujero en la superficie NURBS original.

BIBLIOGRAFÍA

[1] Gerald Farin. Curves and surfaces for CADG. ACADEMIC PRESS, INC., Third Edition, 1973c1993

[2] Foley and van Dam and Feiner and Hughes. Computer Graphics. Principles and practice. Addison-Wesley, second Edition, 1993

[3] R. Löhner and P. Parikh. Generation of three-dimensional unstructured grids by the advancing front method. Int. j. numer. methods fluids. 8, 1988

[4] J. Peraire and J. Peiró and L. Formaggia, K. Morgan and O.C. Zienkiewicz. Finite element Euler computations in three dimensions. Int. j. Numer. Methods Eng. 26, 1988

[5] Peiró J.. A finite element procedure for the solution of the Euler equations on unstructured meshes, 1989

[6] P.L. George and H. Borouchaki. Delaunay Triangulation and Meshing. HERMES, first Edition, 1998

[7] Buscaglia, Gustavo C. and Enzo A. Dari. Anisotropic Mesh Optimization and its Application in Adaptivity. Int. j. numer. methods in eng. 40, 1997

[8] Peter Möller and Peter Hansbo. On advancing front mesh generation in three dimensions. Int. Journal for Num. Meth. in Eng. 38, 1995

[9] Ralf Kreiner. Generation of unstructured Finite Element meshes, 1995

[10] Jas Frykestig. Advancing front mesh generation techniques with application to the Finite Element method, 1994

[11] Rassineux, A. and J-M. Savignat and O. Stab and P. Villon. Surface Remeshing by Local Hermite Diffuse Interpolation. Universit, de Technologie de CompiSgn, 1999

[12] PDA Engineering. PATRAN Plus Users Manual. PATRAN division, Costa Mesa, California, USA, 1991

[13] Altair HyperMesh. http://www.altair.com, 2000

[14] Femsys Femgen/Femview. http://www.femsys.co.uk, 2000

[15] Femap a division of SDRC. http://www.femap.com, 2000

[16] ICEM CFD Engineering. http://www.icemcfd.com, 2000

[17] Ramon Ribó and others. GiD reference manual. CIMNE, 1998

[18] Eugenio Oñate. Cálculo de estructuras por el método de los elementos finitos. CIMNE, 1995

[19] Hoschek, Josef and Dieter Lasser. Fundamentals of computer aided geometric design. Addison-Wesley, 1993

[20] R. Löhner. Regridding surface triangulations. Journal of Computational physics 126, 1996

[21] P. L. George. Automatic Mesh Generation. Application to Finite Element Methods. Wiley, 1991

[22] J. Bonet and J. Peraire. An alternating digital tree (ADT) algorithm for 3D geometric search and itersection problems. Int. Journal for Num. Meth. in Eng. 31, 1991

[23] R. Löhner. Automatic Unstructured Grid Generators. Finite Elements in Analysis and Design 25, 1997

[24] R. Löhner. Some useful data structures for the generation of unstructured grids. Commun. appl. numer. methods 4, 1988

[25] M. Chiumenti and C. Agelet de Saracibar and M. Cervera. Constitutive Modeling and Thermo-Mechanical Phase-Change Systems. CIMNE, 1999

[26] J. Rojek and G. Garcia Garino and E. Oñate. Advanced Finite Element Models for Analysis of Industrial Sheet Forming Processes. CIMNE, 1995

[27] R. Codina. A Finite Element Formulation for Viscous Incompressible Flows. CIMNE, 1993

[28] E. Oñate and S. Oller and S. Botello and J. Miquel Canet. Métodos Avanzados de Cálculo de Estructuras de Materiales Compuestos. CIMNE, 1991

[29] J. Miquel Canet and E. Oñate and C. García Garino and S. Botello and F. Flores and J. Rojek. Análisis de Problemas de Choque e Impacto entre Sólidos Deformables por el Método de los Elementos Finitos.. CIMNE, 1994

[30] F. Zarate and E. Oñate. CALTEP: Programa para el cálculo transitorio de la ecuación de Poisson. CIMNE, 1993

[31] J. Mora, J. Miquel Canet y E. Oñate. Cálculo de Campos Electrostáticos por el Método de Elementos Finitos. CIMNE, 1998

[32] Julio García. A finite element method for hydrodynamic analysis of naval structures, 1999

[33] P. Lancaster and K. Salkauskas. Curve and Surface Fitting. Academic press, first Edition, 1988

[34] J. Peraire and M. Vahdati and K. Morgan and O.C. Zienkiewicz. Adaptive remeshing for compressible flow computations. J. Comp. Phys. 72, 1987

[35] A. Shostko and R. Löhner. Three-Dimensional Parallel Unstructured Grid Generation. Int.J.Num.Meth.Eng. 38, 1995

[36] O. C. Zienkiewicz and R. Taylor. El Método de los Elementos Finitos. Formulación Básica y Problemas Lineales, Volume I. McGraw-Hill, 4${\textstyle ^{a}}$ Edition, 1993

[37] O. C. Zienkiewicz and R. Taylor. El Método de los Elementos Finitos. Mecánica de Sólidos y Fluidos. Dinámica y no Linealidad, Volume II. McGraw-Hill, 4${\textstyle ^{a}}$ Edition, 1994

[38] D.N.Knuth. The art of computer programming, Vol. 3. Addison-Wesley, 1973

[39] NAFEMS. Integrate CAD and Analysis, 1996

[40] J M Le Gouez and M A Boheas and C J Miles. Data requirements for analysis modelling. NAFEMS, 1995

[41] Adi Levin. Interpolating nets of curves by smooth subdivision surfaces. Computer Graphics Proceedings (SIGGRAPH), 1999