Semi-explicit residual goal-oriented estimate for linear mechanics

==Resumen==

El objetivo de este trabajo es definir un método semiexplícito para estimar el error en cantidades de interés derivado de la resolución por elementos finitos de un problema de elasticidad lineal. En este proceso, se identifican dos etapas: una estimación implícita del error del problema adjunto y otra, explícita, asociada al problema directo (primal). La parte implícita de la estimación requiere dos fases, cada una de las cuales consiste en hacer proyecciones sobre espacios de funciones burbuja. Estas proyecciones sobre espacios burbuja comportan operaciones locales y, por tanto, son de bajo coste. Las dos fases se encargan de estimar el error interior en los elementos y la contribución de las aristas (asociada a los saltos de tensiones normales). El carácter implícito de esta parte de la estimación se aplica a la resolución de los problemas locales en espacios de dimensión pequeña. Se analiza también la posibilidad de escoger espacios de funciones burbuja locales de una sola dimensión (fijando la forma de la burbuja y determinando únicamente el coeficiente escalar que la multiplica), lo cual, en la práctica, convierte esta fase en explícita. La parte explícita de la estimación inyecta en el residuo débil primal la aproximación del error adjunto estimada en la primera fase. Esta estrategia es similar a la empleada en el método DWR, pero utilizando una formulación débil del residuo e inyectando un error dual en vez de utilizar una solución posprocesada. La metodología propuesta se valida con su aplicación a un ejemplo numérico.

==Abstract==

We aim at defining a semi-explicit approach to estimate the error in quantities of interest associated with the Finite Element solution of a linear elasticity problem. The advocated procedure is split in two parts, an implicit error estimate for the adjoint problem and an explicit estimate assessing the error in the direct (primal) problem. The implicit part of the estimate (on the adjoint problem) embraces two phases, each consisting in projecting the error on “bubble” functional spaces. The projections are low-cost operations due to their local nature. The two phases account the error in the interior of the elements and the contribution of the elementary edges (associated with the tractions jump). The implicit character of this part is provided by the solution of the local problems in low-dimensional functional spaces. We also analyze the particular case of selecting one-dimensional functional spaces (setting the shape of the bubble and computing a scalar coefficient), which, in practice, make this part of the process explicit. The explicit part of the estimate consists in injecting in the weak primal residual the approximation of the adjoint error obtained in the first phase. This approach is similar to the DWR method but using a weak formulation of the residual and injecting an implicitly estimated adjoint error rather than a post-processed solution. The proposed methodology is validated in a numerical example

==Palabras clave==

Estimación del error ; Semiexplícita ; Funciones burbuja ; Estimadores residuales

==Keywords==

Error estimation ; Semi-explicit ; Bubble functional spaces ; Residual estimates

==1. Introducción==

Cuando se resuelve un problema con el método de los elementos finitos, la solución numérica es una aproximación de la solución exacta, que es desconocida. Existen diferentes métodos para estimar el error cometido (v. [[#bib0105|[1]]]  y sus referencias). Este trabajo se centra en el estudio y el diseño de estimadores a posteriori (frente a los estimadores a priori, que determinan las propiedades de convergencia del método pero no dan un valor preciso del error que se ha cometido realmente) en cantidades de interés (por oposición a los estimadores en norma energética que son los más clásicos). Pensando en extender su dominio de aplicación a problemas de dinámica estructural, interesará utilizar estimadores cuyo coste computacional sea bajo. Por ello, este estudio renuncia, de entrada, a obtener estimadores de error que proporcionen cotas superiores o inferiores del error (v. [[#bib0110|[2]]] , [[#bib0115|[3]]]  y [[#bib0120|[4]]] ) y se conforma con lograr aproximaciones del error razonables. Se busca obtener un estimador que sea rápido y económico en términos de coste computacional, y que pueda ser utilizado en problemas que requieran estimar el error en múltiples ocasiones sobre un mismo dominio, como ocurre en los procesos adaptativos [[#bib0125|[5]]]  o en los problemas transitorios [[#bib0130|[6]]] .      

Para estimar el error en cantidades de interés definidas por el usuario, se introduce un problema auxiliar (adjunto o dual) que permite representar la cantidad de interés en términos de medidas energéticas del problema original (directo o primal) y el adjunto [[#bib0135|[7]]] , [[#bib0140|[8]]]  and [[#bib0145|[9]]] . Las cantidades de interés se describen con funcionales lineales que se utilizan como término fuente (cargas aplicadas) en el problema adjunto. En este trabajo, se presentan representaciones residuales, las cuales utilizan la solución del problema original o primal y de un problema auxiliar o dual para estimar el error. Se muestran las representaciones exactas del error en la cantidad de interés y se proponen estimaciones utilizando varias funciones de tipo burbuja, las cuales se centran en el análisis del error en los desplazamientos y no en las tensiones.      

Los estimadores de tipo ''recovery''  se basan en el posprocesamiento de la solución numérica. Se reconstruye una solución incrementando su regularidad y mejorando su precisión, lo cual sirve para estimar el error cometido en la solución numérica original. La ''super-convergent patch recovery'' , introducida en  [[#bib0150|[10]]] , es una de las opciones más populares y utilizadas para estimar el error en la norma energética y también en cantidades de interés. Aunque, en principio, no proporciona cotas del error, se han introducido variantes que sí las proporcionan (v. [[#bib0155|[11]]]  y [[#bib0160|[12]]] ). Se trata de una familia de estimadores eficiente y robusta. Sin embargo, su extensión a los problemas de dinámica es más compleja, puesto que, en sus versiones originales, no contempla reconstruir campos de desplazamientos o velocidades, únicamente de tensiones.      

En este artículo, se estudia el posible uso de los estimadores residuales en este contexto, con vistas a su generalización al problema dinámico. Se opta por una estrategia residual semiexplícita, en la cual se espera que el coste sea moderado. En la perspectiva de la dinámica estructural, renunciar a las propiedades de acotación que proporcionan los estimadores implícitos no es grave (en todo caso, obtener cotas es muy costoso y estas no son muy precisas cuando el amortiguamiento es pequeño [[#bib0165|[13]]] ). El carácter semiexplícito se asocia a una aproximación implícita del error en el problema adjunto y a una estimación explícita del error en la cantidad de interés a partir de una representación del error que consiste en evaluar el residuo débil primal (explícito) en la estimación del adjunto (o dual). Esta naturaleza semiexplícita hace que el método sea poco costoso, y se prevé extenderlo a problemas de dinámica estructural en que el coste computacional es un cuello de botella. En la resolución de la parte implícita (estimación del problema adjunto), y también para el cálculo del estimador explícito, es necesario introducir una malla refinada en que se calculan las aproximaciones del error y de los residuos. Se consideran dos estrategias distintas para construir esta malla refinada: una con una submalla estructurada y otra, no estructurada. Existen varias formulaciones de estimadores residuales implícitos, todas ellas basadas en resolver problemas locales con diferentes estrategias para descomponer el dominio e imponer las condiciones de contorno correspondientes (v. [[#bib0170|[14]]]  y [[#bib0175|[15]]] ). En este trabajo, se opta por construir la estimación con funciones positivas de soporte compacto, denominadas ''funciones de tipo burbuja''  (v.  [[#bib0170|[14]]]  y [[#bib0180|[16]]] ), y analizar el error relativo y el índice de efectividad correspondiente.      

La estructura del resto del artículo es la siguiente: en la sección 2, se presenta el problema a resolver, así como la forma de medir los errores de aproximación. En la sección 3, se introducen la cantidad de interés y el problema auxiliar utilizado en las representaciones del error que se dan en la sección siguiente. En la sección 4, se analiza un nuevo término que aparece en las representaciones del error en cantidades de interés cuando se utilizan submallas distorsionadas. En la sección 5, se presenta una estimación residual implícita del error [[#bib0120|[4]]] , para la cual se requiere resolver un problema sobre cada elemento del dominio y un problema por cada arista. Para reducir el coste computacional de este procedimiento, en la sección 6 se presentan dos estimaciones del error utilizando funciones de tipo burbuja. Los resultados numéricos se presentan en la sección 7 y las conclusiones, en la sección 8.      

==2. Planteamiento del problema==

Sea <math display="inline">\Omega \subset R^2</math>  un dominio abierto y acotado con frontera continua y función lipschitziana <math display="inline">\partial \Omega </math> , compuesta por una parte del tipo Dirichlet y otra del tipo Neumann, es decir, <math display="inline">\partial \Omega =\overline{\Gamma _D}\cup \overline{\Gamma _N}</math> , con <math display="inline">\overline{\Gamma _D}\cap \overline{\Gamma _N}=\varnothing </math> . El dominio está ocupado por un material elástico lineal, que está sometido a la acción de una fuerza '''f'''  en <math display="inline">[L^2(\Omega )]^2</math>  y a tensiones de borde '''h'''  en <math display="inline">\left[L^2\left(\Gamma _N\right)\right]^2</math> . Los desplazamientos <math display="inline">\mbox{u}</math>  y las tensiones <math display="inline">\sigma =\sigma \left(u\right)</math>  en el cuerpo elástico satisfacen la condición de equilibrio y las condiciones de contorno siguientes:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>-\nabla \cdot \sigma \left(\mbox{u}\right)=\mbox{f}\mbox{ }\mbox{sobre}\mbox{ }\Omega </math>
|}
| style="text-align: right;" | ( 1.1)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{u}=0\mbox{ }\mbox{en}\mbox{ }\Gamma _D</math>
|}
| style="text-align: right;" | ( 1.2)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{n}\cdot \sigma (\mbox{u})=\mbox{h}\mbox{ }\mbox{en}\mbox{ }\Gamma _N</math>
|}
| style="text-align: right;" | ( 1.3)
|}

donde '''n'''  denota la normal unitaria exterior a <math display="inline">\Gamma _N</math> . El tensor de tensiones <math display="inline">\sigma (\mbox{u})</math>  está relacionado con el tensor de deformaciones

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\epsilon \left(\mbox{u}\right)=\frac{1}{2}\left(\nabla \mbox{u}+\nabla ^T\mbox{u}\right)\equiv \left\{\epsilon _x\epsilon _y\gamma _{xy}\right\}^T</math>
|}
| style="text-align: right;" | 
|}

por medio de la ley de Hooke

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\sigma \left(\mbox{u}\right)=\mbox{C}:\epsilon \left(\mbox{u}\right)\equiv \left\{\sigma _x\sigma _y\tau _{xy}\right\}^T\mbox{,}</math>
|}
| style="text-align: right;" | 
|}

donde '''C'''  es el tensor de constantes elásticas. Con el método de los residuos ponderados de Galerkin, se obtiene la forma variacional del problema, que consiste en encontrar una función '''u'''  en el espacio V  tal que

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a\left(\mbox{u},\mbox{v}\right)=l\left(\mbox{v}\right),\forall \mbox{v}\in \mbox{V,}</math>
|}
| style="text-align: right;" | ( 2)
|}

donde <math display="inline">a(\mbox{u},\mbox{v})=\int_{\Omega }\sigma \left(\mbox{u}\right):\epsilon \left(\mbox{v}\right)d\Omega </math> , <math display="inline">l\left(\mbox{v}\right)=\int_{\Omega }\mbox{f}\cdot \mbox{v}d\Omega +\int_{\Gamma _N}\mbox{h}\cdot \mbox{v}d\Gamma </math>  y <math display="inline">\mbox{V}=\left\{\mbox{w}\in \left(H^1\left(\Omega \right)\right)^2:\mbox{w}=0\mbox{ }en\mbox{ }\Gamma _D\right\}</math> , con <math display="inline">H^1(\Omega )</math> , que denota un espacio de Sobolev. La norma en energía se define por

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left|\left|\mbox{v}\right|\right|_E:=\sqrt{a\left(\mbox{v},\mbox{v}\right)}\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

Para aproximar la solución '''u''' , se construye un espacio de elementos finitos V<sub>H</sub> , inducido por una partición del dominio <math display="inline">\Omega </math>  en elementos <math display="inline">\Omega _k,k=1,\cdots ,n_{elem}\mbox{,}</math>  tales que <math display="inline">\overline{\Omega }=\underset{k}{\cup }\overline{\Omega }_k</math> , con <math display="inline">\Omega _k\cap \Omega _j=\varnothing </math> , para <math display="inline">k\not =j</math> . La solución de elementos finitos <math display="inline">\mbox{u}_H</math>  pertenece al espacio V<sub>H</sub>  y satisface la ecuación

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a\left(\mbox{u}_H,\mbox{v}_H\right)=L\left(\mbox{v}_H\right),\forall \mbox{v}_H\in V_H\mbox{.}</math>
|}
| style="text-align: right;" | ( 3)
|}

El error numérico asociado a <math display="inline">\mbox{u}_H</math> , se denota por <math display="inline">\mbox{e}</math>  y está definido por <math display="inline">\mbox{e}=\mbox{u}-\mbox{u}_H</math> . Este error pertenece al espacio V  cuando <math display="inline">u_H</math>  es interpolado en V  y satisface la ecuación del error

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a\left(\mbox{e},\mbox{v}\right)=L\left(\mbox{v}\right)-a\left(\mbox{u}_H,\mbox{v}\right)=:R^P\left(\mbox{v}\right),\forall \mbox{v}\in V\mbox{,}</math>
|}
| style="text-align: right;" | ( 4)
|}

donde <math display="inline">R^P</math>  se denomina ''residuo débil  ''  asociado a <math display="inline">\mbox{u}_H</math>  o ''residuo primal'' . Cuando no se dispone de una solución exacta del problema, se utiliza una malla más fina, denominada ''malla de referencia  '' , para obtener una solución de referencia, <math display="inline">\mbox{u}_{ref}</math> . Con esta solución, el error se define por <math display="inline">\mbox{e}_{ref}=\mbox{u}_{ref}-\mbox{u}_H</math> , donde <math display="inline">\mbox{u}_H</math>  es interpolada sobre la malla de referencia.      

==3. Cantidad de interés y problema dual==

Sea <math display="inline">J</math>  un funcional lineal. Se quiere calcular el error en la evaluación de <math display="inline">J(\mbox{u}_H)</math> . Para ello, se introduce un problema auxiliar, que consiste en encontrar una función <math display="inline">\mbox{u}^\ast </math>  en el espacio V , tal que

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a(\mbox{v},\mbox{u}^\ast )=J(\mbox{v}),\forall \mbox{v}\in V\mbox{.}</math>
|}
| style="text-align: right;" | ( 5)
|}

El error numérico asociado a una solución de elementos finitos <math display="inline">\mbox{u}_H^\ast </math>  de este problema dual, calculada sobre V<sub>H</sub> , se denota por <math display="inline">\epsilon </math>  y está definido por <math display="inline">\epsilon =\mbox{u}^\ast -\mbox{u}_H^\ast </math> . Este error pertenece al espacio V  cuando <math display="inline">\mbox{u}_H^\ast </math>  es interpolado en V  y satisface la ecuación del error

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a\left(\epsilon ,\mbox{v}\right)=J\left(\mbox{v}\right)-a\left(\mbox{v},\mbox{u}_H^\ast \right)=:R^D\left(\mbox{v}\right),\forall \mbox{v}\in V\mbox{,}</math>
|}
| style="text-align: right;" | ( 6)
|}

donde <math display="inline">R^D</math>  se denomina ''residuo débil  ''  asociado a <math display="inline">\mbox{u}_H^\ast </math>  o ''residuo dual  '' . Una de las cantidades de interés analizadas en este trabajo es un valor puntual de la solución. Se quiere utilizar la solución numérica del problema para obtener un valor suficientemente preciso de la solución en un punto <math display="inline">\mbox{x}_0\in \Omega </math>  dado. Para lograrlo, se utiliza el funcional

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>J\left(\mbox{u}\right)=\int_{\Omega }\mbox{q}\cdot \mbox{u}W_r\left(\mbox{x}-\mbox{x}_0\right)d\Omega </math>
|}
| style="text-align: right;" | ( 7)
|}

donde <math display="inline">\mbox{q}=\left\{q_xq_y\right\}^T</math>  representa un vector que indica la dirección y el sentido de una fuerza en el problema dual. La magnitud <math display="inline">J(\mbox{u})(\simeq q_xu_x\left(\mbox{x}_0\right)+q_yu_y(\mbox{x}_0)</math> ) representa el promedio de los desplazamientos '''u'''  en la dirección de <math display="inline">\mbox{q}</math> , en un pequeño entorno de <math display="inline">x_0</math> . La función de suavizado, <math display="inline">W_r</math>[[#bib0185|[17]]]  se escoge del modo siguiente:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>W_r(x)=\left\{\begin{array}{cc}
c\mbox{ }exp\left(\frac{-r^2}{r^2-\left|\left|x\right|\right|^2}\right) & si\left|\left|x\right|\right|<r\\
0 & en\mbox{ }otro\mbox{ }caso
\end{array}\right.</math>
|}
| style="text-align: right;" | 
|}

con <math display="inline">r>0</math> , y la constante <math display="inline">c=c(r,x_0)</math>  se escoge de modo que <math display="inline">\int_{\Omega }W_r\left(x-x_0\right)d\Omega =1.</math>

==4. Representaciones residuales del error==

Si se dispone de la solución exacta del problema, el error en la cantidad de interés descrita a través del funcional <math display="inline">J</math>  se puede representar de tres formas distintas: utilizando el producto escalar de la energía de deformación, a través del residuo primal y a través del residuo dual. Pero, si no se dispone de la solución exacta, se puede encontrar una solución de referencia, calculada sobre una malla suficientemente fina, que sea una buena aproximación de la solución y, en consecuencia, del error exacto. En este caso, si con la malla fina no se genera un subespacio funcional que contenga el subespacio asociado a la malla gruesa, aparecerá un nuevo término en las representaciones del error en cantidades de interés.      

===4.1. Representaciones exactas del error===

Si <math display="inline">\mbox{u}_H</math>  y <math display="inline">\mbox{u}_H^\ast </math>  son aproximaciones en <math display="inline">\mbox{V}_\mbox{H}</math>  de los problemas primal y dual dados en las ecuaciones (2) y (5), el error en la cantidad de interés asociado a estas soluciones se puede representar de las tres formas siguientes [[#bib0135|[7]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>J\left(\mbox{e}\right)=a(\mbox{e},\epsilon )</math>
|}
| style="text-align: right;" | ( 8)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>=R^P\left(\epsilon \right)</math>
|}
| style="text-align: right;" | ( 9)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>=R^D(\mbox{e})</math>
|}
| style="text-align: right;" | ( 10)
|}

donde <math display="inline">a(\ast ,\ast )</math>  es el producto escalar ''energía de deformación  '' ; <math display="inline">R^P\left(\ast \right)</math> , el residuo primal, y <math display="inline">R^D(\ast )</math> , el residuo dual, dados en las ecuaciones (2), (4) y (6). Para obtener estas representaciones, se parte de la validez de la ortogonalidad de Galerkin,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a\left(\mbox{e},\mbox{u}_H^\ast \right)=R^P\left(\mbox{u}_H^\ast \right)=0\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

Como <math display="inline">e\in V</math> , escogiendo <math display="inline">\mbox{v}=\mbox{e}</math>  en la ecuación (5),

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>J\left(\mbox{e}\right)=a\left(\mbox{e},\mbox{u}^\ast \right)\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

Teniendo en cuenta que <math display="inline">\mbox{u}^\ast =\mbox{u}_H^\ast +\epsilon </math>  y la ortogonalidad de Galerkin, se sigue que

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>J\left(\mbox{e}\right)=a\left(\mbox{e},\mbox{u}_H^\ast \right)+a\left(\mbox{e},\epsilon \right)=R^P\left(\mbox{u}_H^\ast \right)+a\left(\mbox{e},\epsilon \right)=a\left(\mbox{e},\epsilon \right)\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

Por otro lado, como <math display="inline">\epsilon \in V</math> , se puede utilizar en la ecuación (4) para obtener

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}\displaystyle R^P\left(\epsilon \right)=L\left(\epsilon \right)-a\left(\mbox{u}_H,\epsilon \right)=a\left(\mbox{u},\epsilon \right)-a\left(\mbox{u}_H,\epsilon \right)=a\left(e\right.\\\left.\displaystyle \epsilon \right)=J\left(e\right)-R^P\left(\mbox{u}_H^\ast \right)=J\left(e\right)\mbox{.}\end{array}</math>
|}
| style="text-align: right;" | 
|}

Finalmente, utilizando la ecuación (6),

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R^D\left(e\right)=J\left(e\right)-a\left(e,\mbox{u}_H^\ast \right)=J\left(e\right)-R^P\left(\mbox{u}_H^\ast \right)=J\left(e\right)\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

===4.2. Análisis del término <math display="inline">R^P\left(\mbox{u}_H^\ast \right)</math>===

Las representaciones del error que se muestran en el apartado 4.1 se basan en la propiedad de ortogonalidad de Galerkin. Cuando se calcula una solución aproximada sobre una malla gruesa y esta es interpolada sobre una malla refinada encajada a la anterior, y, además, los espacios funcionales asociados a ambas mallas están encajados, la solución interpolada satisface la propiedad. Sin embargo, si no se cumple la contención de los espacios, el término <math display="inline">R^P\left(\mbox{u}_H^\ast \right)</math>  puede ser no nulo, debido a que, al calcular las integrales involucradas en dicho término, se generan errores, ya que la solución interpolada no puede representarse de forma exacta sobre la malla fina. Este fenómeno se observa cuando se utilizan funciones bilineales típicas para elementos cuadriláteros. En esta sección, se analiza el término <math display="inline">R^P\left(\mbox{u}_H^\ast \right)</math> , calculado sobre una malla distorsionada de elementos cuadriláteros encajada en una malla gruesa. Se quiere ver el efecto de la distorsión de una malla refinada en el cálculo de la integral de una función de forma proyectada sobre dicha malla y su consecuencia sobre el término <math display="inline">R^P\left(\mbox{u}_H^\ast \right)</math> , dado en las representaciones del error en cantidades de interés.      

Sea <math display="inline">N_1\left(\xi ,\eta \right)=0,25\mbox{ }(1-\xi )(1-\eta )</math>  una función de forma típica en el cálculo de elementos finitos, definida sobre el cuadrado <math display="inline">\left[-1,1\right]\times [-1,1]</math>  en el sistema local <math display="inline">\xi \eta </math> . Esta función se puede enviar a <math display="inline">\Omega =\left[-1,1\right]\times \left[-1,1\right]</math>  en un sistema global <math display="inline">xy\mbox{,}</math>  teniendo en este la representación <math display="inline">N_1\left(x,y\right)=0,25\left(1-x\right)\left(1-y\right)</math> . Calculando, se obtiene que <math display="inline">\int_{\Omega }N_1d\Omega =0,25</math> .      

Si se escoge una partición de <math display="inline">\Omega </math>  que tenga cuatro elementos con sus aristas ortogonales, <math display="inline">\Omega =\underset{\mbox{j}=1}{\overset{4}{\cup }}\Omega _j</math> , entonces, como consecuencia de la contención de los subespacios asociados a ambas mallas y de la bilinealidad de las <math display="inline">N_i</math> , se sigue que

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\int_{\Omega }N_1\mbox{d}\Omega =\sum_{\mbox{j}=1}^4\int_{\Omega _\mbox{j}}N_1\mbox{d}\Omega =\frac{1}{4}.</math>
|}
| style="text-align: right;" | 
|}

Por el contrario, si se escoge una partición del dominio formada por elementos distorsionados, como la que se muestra en la [[#fig0005|figura 1]] , el valor de la integral será distinto de <math display="inline">\frac{1}{4}</math> . En la [[#fig0010|figura 2]] , se pueden observar los contornos de las superficies correspondientes a la parte de <math display="inline">N_1</math>  definida sobre los elementos de la malla distorsionada. La traza de la gráfica de la función de forma <math display="inline">N_1</math>  definida sobre el contorno distorsionado está formada por arcos de curvas (cuadráticas), mientras que la traza correspondiente a la proyección está formada por segmentos de recta. Luego, al calcular la integral de <math display="inline">N_1</math>  utilizando dicha aproximación, se generan errores de proyección. Observando la [[#fig0010|figura 2]] , se puede inferir que, al utilizar una malla refinada con elementos distorsionados, aunque sea encajada, el espacio asociado a la malla distorsionada no va a contener al espacio asociado a la malla grosera, pues no se puede representar <math display="inline">N_1</math>  de forma exacta sobre la malla fina. Esto, a su vez, hace que, en la representación del error en la cantidad de interés, se genere un nuevo error, que es captado exactamente con el término <math display="inline">R^P\left(\mbox{u}_H^\ast \right)</math> , como se evidencia en los resultados numéricos.

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_530383077-1-s2.0-S0213131515000504-gr1.jpg|center|349px|Malla con elementos distorsionados.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

Malla con elementos distorsionados.

</span>
|}

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_530383077-1-s2.0-S0213131515000504-gr2.jpg|center|344px|Diferencias entre los contornos de la gráfica de la función y su interpolación ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

Diferencias entre los contornos de la gráfica de la función y su interpolación sobre los elementos de una malla distorsionada.

</span>
|}

Cabe destacar que esta situación no se presenta cuando la malla refinada está formada por elementos no distorsionados. No obstante, en procesos de adaptatividad, utilizando elementos cuadriláteros con funciones de forma bilineales, se pueden generar submallas con elementos distorsionados, lo cual justifica el análisis anterior.

===4.3. Nuevas representaciones del error===

Al proyectar una solución sobre una malla refinada que tenga elementos distorsionados, como la que se puede observar en la [[#fig0005|figura 1]] , se generan errores en el cálculo del término <math display="inline">R^P\left(\mbox{u}_H^\ast \right)</math> . Estos errores inciden en la representación del error en cantidades de interés, lo cual se expresa en las representaciones siguientes:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>J(e)=a(e,\epsilon )+R^P(\mbox{u}_H^\ast )</math>
|}
| style="text-align: right;" | ( 11)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>=R^P\left(\epsilon \right)+R^P\left(\mbox{u}_H^\ast \right)=:Rep_2</math>
|}
| style="text-align: right;" | ( 12)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>=R^D\left(e\right)+R^P\left(\mbox{u}_H^\ast \right)=:Rep_1\mbox{.}</math>
|}
| style="text-align: right;" | ( 13)
|}

Aquí, las soluciones aproximadas al problema, así como la de <math display="inline">\mbox{u}_H^\ast </math> , son interpoladas sobre la malla de referencia utilizada para calcular los errores.      

Para verificar la primera representación, en la ecuación (5) se toma <math display="inline">\mbox{u}^\ast =\mbox{u}_H^\ast +\epsilon </math> , de donde <math display="inline">a\left(\mbox{v},\mbox{u}_H^\ast +\epsilon \right)=J\left(\mbox{v}\right)</math> . Luego, por la linealidad de <math display="inline">a\left(\ast ,\ast \right)</math> , <math display="inline">J(\mbox{v})=a(\mbox{v},\mbox{u}_H^\ast )+a(\mbox{v},\epsilon )</math> . Escogiendo <math display="inline">\mbox{v}=e</math>  en la última ecuación, se tiene <math display="inline">J(e)=a(e,\mbox{u}_H^\ast )+a(e,\epsilon )</math> , pero

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}\displaystyle a\left(e,\mbox{u}_H^\ast \right)=a\left(\mbox{u}-\mbox{u}_H,\mbox{u}_H^\ast \right)=L\left(\mbox{u}_H^\ast \right)-a\left(\mbox{u}_H\right.\\\left.\displaystyle \mbox{u}_H^\ast \right)=R^P\left(\mbox{u}_H^\ast \right)\mbox{.}\end{array}</math>
|}
| style="text-align: right;" | 
|}

Por tanto,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>J\left(e\right)=a\left(e,\epsilon \right)+R^P(\mbox{u}_H^\ast )\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

Para verificar la segunda representación, se escoge <math display="inline">\mbox{v}=\epsilon </math>'''e''' n la ecuación (4):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}\displaystyle R^P\left(\epsilon \right)=L\left(\epsilon \right)-a\left(\mbox{u}_H,\epsilon \right)=a\left(\mbox{u},\epsilon \right)-a\left(\mbox{u}_H,\epsilon \right)=a\left(e\right.\\\left.\displaystyle \epsilon \right)=J\left(e\right)-R^P\left(\mbox{u}_H^\ast \right)\mbox{.}\end{array}</math>
|}
| style="text-align: right;" | 
|}

Por tanto,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>J\left(e\right)=R^P\left(\epsilon \right)+R^P\left(\mbox{u}_H^\ast \right)\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

Para verificar la tercera representación, se parte de la ecuación (6): <math display="inline">R^D\left(\mbox{v}\right)=J\left(\mbox{v}\right)-a\left(\mbox{v},\mbox{u}_H^\ast \right)</math> , donde se escoge <math display="inline">\mbox{v}=e</math> , para obtener

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R^D\left(e\right)=J\left(e\right)-a\left(e,\mbox{u}_H^\ast \right)=J\left(e\right)-R^P\left(\mbox{u}_H^\ast \right)\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

Por tanto,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>J\left(e\right)=R^D\left(e\right)+R^P\left(\mbox{u}_H^\ast \right)\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

Las representaciones residuales del error presentadas anteriormente permiten utilizar cualquier estimador del error en los desplazamientos para calcular aproximadamente el error en la cantidad de interés.

==5. Estimación implícita del error==

Una forma de calcular el error de una solución de elementos finitos es resolver problemas locales para aproximar las dos contribuciones, la debida al error en el interior de los elementos y la debida al error sobre las aristas. La suma de las soluciones de estos problemas locales se denota por <math display="inline">\overline{e}=\overline{e}^{int}+\overline{e}^{edg}\approx e</math> .      

Para calcular aproximadamente los errores interiores, se escribe <math display="inline">\overline{e}^{int}\approx \sum_{k=1}^{n_{elem}}\overline{e}_k^{int}</math> , donde <math display="inline">\overline{e}_k^{int}</math>  representa la contribución del elemento <math display="inline">\Omega _k</math>  al error interior. Luego, utilizando la ecuación del error (4) sobre cada uno de los elementos de la malla, se resuelve un problema con condiciones de contorno de Dirichlet homogéneas, tomando como término fuente la restricción del residuo débil al elemento <math display="inline">\Omega _k</math> . El problema local consiste en resolver, sobre una malla refinada como la que se muestra en la [[#fig0015|figura 3]] , la ecuación

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a_k\left(\overline{e}_k^{int},v\right)=R_k^P\left(\mbox{v}\right)\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_530383077-1-s2.0-S0213131515000504-gr3.jpg|center|357px|Malla gruesa, malla refinada y dominio del problema local para los errores ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 3.

Malla gruesa, malla refinada y dominio del problema local para los errores interiores.

</span>
|}

Los operadores <math display="inline">a_k\left(\ast ,\ast \right)</math>  y <math display="inline">R_k^P\left(\ast \right)</math>  representan las restricciones de los operadores correspondientes al elemento <math display="inline">\Omega _k</math> . Para calcular aproximadamente el error debido a la contribución sobre las aristas, se escribe <math display="inline">\overline{e}^{edg}\approx \sum_{l=1}^{n_{edg}}\overline{e}_l^{edg}</math> , donde <math display="inline">\overline{e}_l^{edg}</math>  representa la contribución de la arista <math display="inline">l</math>  al error y <math display="inline">n_{edg}</math>  es el número de aristas. Para cada una de las aristas, se escoge un subdominio <math display="inline">\Lambda _l\subset \Omega </math>  que la contenga, como el representado en la [[#fig0020|figura 4]] . Este subdominio estará contenido en dos elementos si la arista es interior y en un solo elemento si es una arista de la frontera de <math display="inline">\Omega </math> . El problema consiste en resolver, sobre una malla refinada como la que se puede observar en la [[#fig0020|figura 4]] , la ecuación:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a_l\left(\overline{e}_l^{edg},\mbox{v}\right)=R_l^P\left(\mbox{v}\right)-a_l\left(\overline{e}^{int},\mbox{v}\right)\mbox{,}</math>
|}
| style="text-align: right;" | 
|}

tomando condiciones de Dirichlet homogéneas sobre el contorno. El error sobre arista encontrado como la solución a dicho problema se ortogonaliza de los errores interiores mediante la fórmula:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}\displaystyle \overline{e}_l^{edg}=\overline{e}_l^{edg}-\frac{a\left(\overline{e}_l^{edg},\overline{e}_k^{int}\right)}{a\left(\overline{e}_k^{int},\overline{e}_k^{int}\right)}e_k^{int}-\\\displaystyle -\frac{a\left(\overline{e}_l^{edg},\overline{e}_{k^'}^{int}\right)}{a\left(\overline{e}_{k^'}^{int},\overline{e}_{k^'}^{int}\right)}e_{k^'}^{int}\mbox{,}\end{array}</math>
|}
| style="text-align: right;" | 
|}

donde los subíndices <math display="inline">k</math>  y <math display="inline">k'</math>  corresponden a los índices de los elementos <math display="inline">\Omega _k</math>  y <math display="inline">\Omega _{k'}</math>  que intersecan el subdominio <math display="inline">\Lambda _l</math> . Una vez calculadas las contribuciones sobre el interior y sobre las aristas, el error estimado se escribe como

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overline{e}:=\sum_{k=1}^{n_{elem}}\overline{e}_k^{int}+\sum_{l=1}^{n_{edg}}\overline{e}_l^{edg}\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_530383077-1-s2.0-S0213131515000504-gr4.jpg|center|357px|Malla gruesa, malla refinada y dominio del problema local para los errores sobre ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 4.

Malla gruesa, malla refinada y dominio del problema local para los errores sobre aristas.

</span>
|}

Cuando se ha estimado el error implícitamente, se utiliza el residuo dual para obtener una estimación del error en la cantidad de interés, según la representación dada en la ecuación (10):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}\displaystyle E_1:=R^D\left(\overline{e}\right)=\sum _{k=1}^{n_elem}R^D\left(\overline{e}_k^{int}\right)+\\\displaystyle +\sum _{l=1}^{n_edg}R^D(\overline{e}_l^{edg})\mbox{.}\end{array}</math>
|}
| style="text-align: right;" | ( 14)
|}

Para tener un mapa del error estimado sobre todo el dominio, la distribución se hace por elementos. El error interior <math display="inline">R^D\left(\overline{\mbox{e}}_k^{int}\right)</math>  es asignado al elemento <math display="inline">\Omega _k</math> . El error sobre arista <math display="inline">R^D(\overline{e}_l^{edg})</math>  es asignado a uno o dos elementos; Si la arista es interior, la contribución se distribuye de forma equitativa entre los dos elementos que intersecan el subdominio <math display="inline">\Lambda _l</math> . Si la arista es frontera, la contribución es asignada completamente al elemento que contiene el subdominio.      

Procediendo de forma análoga a la descrita anteriormente, se obtiene otra estimación del error en la cantidad de interés, pero utilizando el error del problema auxiliar:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{c}\overline{\epsilon }=\overline{\epsilon }^{int}+\overline{\epsilon }^{edg}\approx \epsilon \\ E_2:=R^P\left(\overline{\epsilon }\right):=\sum _{k=1}^{n_elem}R^P\left(\overline{\epsilon }_k^{int}\right)+\\\displaystyle +\sum _{l=1}^{n_edg}R^P\left(\overline{\epsilon }_l^{edg}\right). \end{array}</math>
|}
| style="text-align: right;" | ( 15)
|}

Para escribir cada una de las dos representaciones anteriores, se recuperan los errores estimados elemento a elemento y se calculan los residuos correspondientes sobre el interior y sobre las aristas.

==6. Estimación con función burbuja==

Se extiende el uso de las funciones de tipo burbuja a problemas vectoriales, como los presentadas en [[#bib0180|[16]]] , donde se estima el error en cantidades de interés utilizando funciones de tipo burbuja para problemas elípticos, siguiendo las ideas de [[#bib0170|[14]]] . Las estimaciones son de tipo residual dual ponderado, como las presentadas en [[#bib0190|[18]]]  y [[#bib0195|[19]]] . En [[#bib0180|[16]]] , se estima el error en cantidades de interés, y el error en desplazamientos se aproxima como la suma de múltiplos de funciones de tipo burbuja, definidas con soporte compacto sobre todos los elementos de una malla de elementos finitos. Esta suma consta de dos contribuciones al error, una sobre el interior de los elementos y otra sobre las aristas. Luego, el error en la cantidad de interés se estima utilizando el error aproximado y los residuos de los problemas primal y dual según las representaciones allí definidas. En este trabajo, se ensaya con varias funciones de tipo burbuja diferentes, para analizar tanto el comportamiento del error relativo como el de los índices de efectividad en cada caso. Se trabaja con dos tipos de funciones con soporte compacto. Las primeras se escogen con derivada no nula en la dirección normal sobre el contorno de cada elemento y las segundas, con derivada nula en la dirección normal sobre el contorno de cada elemento. El análisis de los resultados obtenidos permitirá determinar cuál de todas las anteriores es la función de tipo burbuja más adecuada para estimar el error.      

===6.1. Estimación interior y sobre aristas con funciones del tipo burbuja===

Para estimar el error en cantidades de interés, se utilizan dos soluciones de elementos finitos. Una correspondiente al problema original o primal (ecuación 2) y otra correspondiente a un problema auxiliar o dual (ecuación 5), que se resuelven sobre una misma malla (aunque se pueden resolver sobre mallas diferentes). Se utilizan las representaciones del error definidas en las ecuaciones (9) y (10), utilizando el residuo de los problemas primal y dual. Los residuos de las representaciones operan sobre una malla refinada o de referencia, que se escoge encajada a la malla gruesa. La estimación de este error con funciones de tipo burbuja se hace en dos etapas. En la primera etapa, se estiman los errores interiores, <math display="inline">e^{int}</math> , correspondientes a las contribuciones al error sobre el interior de los elementos. Esta aproximación se denota por <math display="inline">\tilde{e}^{int}</math> . Para ello, se supone que, sobre cada elemento, el error interior se puede escribir como un múltiplo de una función burbuja, de soporte compacto sobre el elemento:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>e_k^{int}\approx \tilde{e}_k^{int}:=c_k\psi _k=c_k\lbrace \psi _k^x,(1+\nu )\mbox{ }\psi _k^y\rbrace \mbox{,}</math>
|}
| style="text-align: right;" | 
|}

en que la función burbuja se da a priori, sin tener que resolver un problema de Dirichlet como en el caso de la estimación implícita; <math display="inline">\nu </math>  es el radio de Poisson, y las constantes <math display="inline">c_k</math>  se han de determinar. Esta hipótesis permite que el cálculo aproximado del error interior consista en realizar solo productos de matrices con vectores. Aunque se perderá un poco la precisión de la estimación del error, el cálculo será mucho más rápido, sobre todo cuando se tenga un gran número de elementos, o cuando haya que calcular el error en muchas ocasiones. Entonces, a partir de la ecuación

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a_k\left(\tilde{e}_k^{int},\mbox{v}\right)=R_k^P\left(\mbox{v}\right)\mbox{,}</math>
|}
| style="text-align: right;" | 
|}

y escogiendo <math display="inline">\mbox{v}=\psi _k</math> , se sigue que

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a_k\left(c_k\psi _k,\psi _k\right)=R_k^P\left(\psi _k\right)\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

de donde se puede obtener <math display="inline">c_k</math>  y el error interior queda expresado como

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tilde{e}_k^{int}=\frac{R_k^P\left(\psi _k\right)}{a_k\left(\psi _k,\psi _k\right)}\psi _k\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

En la segunda etapa, se estima la contribución al error sobre las aristas, <math display="inline">\mbox{e}^{edg}</math> . Esta aproximación se denota por <math display="inline">\tilde{\mbox{e}}^{edg}</math> . Para ello, sobre cada subdominio <math display="inline">\Lambda _l</math>  se escoge una función burbuja <math display="inline">\chi _l</math>  de soporte compacto y se supone que, sobre dicho subdominio, el error sobre arista se aproxima como un múltiplo de la burbuja. Esto es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>e_l^{edg}\approx \tilde{e}_l^{edg}:=d_l\chi _l=d_l\left\{\chi _l^x,\left(1+\nu \right)\chi _l^y\right\}\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

En la ecuación <math display="inline">a\left(e^{int}+e^{edg},\mbox{v}\right)=R^P\left(\mbox{v}\right)</math> , se utilizan las aproximaciones del error <math display="inline">a\left(\tilde{e}^{int}+\tilde{e}^{edg},\mbox{v}\right)=R^P\left(\mbox{v}\right)</math> . Por la linealidad del producto escalar, se sigue que <math display="inline">a\left(\tilde{e}^{edg},\mbox{v}\right)=R^P\left(\mbox{v}\right)-a\left(\tilde{e}^{int},\mbox{v}\right)</math> . Luego, escogiendo <math display="inline">\mbox{v}=\chi _l</math>  y teniendo en cuenta la ortogonalidad de las funciones burbuja (por ser escogidas con soporte compacto), se tiene que

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a\left(d_l\chi _l,\chi _l\right)=R^P\left(\chi _l\right)-a\left(\tilde{e}^{int},\chi _l\right)\mbox{,}</math>
|}
| style="text-align: right;" | 
|}

de donde se puede obtener <math display="inline">d_l</math>  y, denotando la restricción del residuo débil al subdominio <math display="inline">\Lambda _l</math>  por <math display="inline">R_l^P\left(\ast \right)</math> , se sigue que

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tilde{e}_l^{edg}=\frac{R_l^P\left(\chi _l\right)-a\left(\tilde{e}^{int},\chi _l\right)}{a\left(\chi _l,\chi _l\right)}\chi _l\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

Una vez obtenida la estimación sobre aristas, se ortogonaliza a partir de la estimación interior con el proceso de Gram-Schmidt. La estimación del error primal, denotada por <math display="inline">\tilde{e}</math> , es definida como la suma de las dos contribuciones, es decir, <math display="inline">e\approx \tilde{e}=\tilde{e}^{int}+\tilde{e}^{edg}</math> . Con esta estimación y la representación dada en la ecuación (10), se obtiene la estimación del error en cantidades de interés:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>E_3:=R^D\left(\tilde{e}\right)=R^D(\tilde{e}^{int})+R^D(\tilde{e}^{edg})\mbox{,}</math>
|}
| style="text-align: right;" | ( 16)
|}

donde <math display="inline">R^D\left(\tilde{e}^{int}\right)=\sum_{k=1}^{n_{elem}}\frac{R^P\left(\psi _k\right)R^D\left(\psi _k\right)}{a\left(\psi _k,\psi _k\right)}</math>  y <math display="inline">R^D\left(\tilde{e}^{edg}\right)=\sum_{l=1}^{n_{edg}}\frac{R^P\left(\chi _l\right)-a\left(\chi _l,\tilde{e}^{int}\right)}{a\left(\chi _l,\chi _l\right)}R^D\left(\chi _l\right)</math>  y las componentes de las funciones <math display="inline">\psi _k=\left\{\psi _{k_x},\left(1+\nu \right)\psi _{k_y}\right\}</math> , para todo <math display="inline">k=1,\ldots ,n_{elem}</math>  y <math display="inline">\chi _l=\left\{\chi _{l_x},\left(1+\nu \right)\chi _{l_y}\right\}\mbox{,}</math>  para todo <math display="inline">l=1,\ldots ,n_{edg}</math> , son del tipo burbuja, definidas sobre el interior y sobre las aristas de los elementos, respectivamente.      

Para estimar el error del problema dual, se procede de manera análoga a como se ha hecho anteriormente. En la primera etapa, se estiman los errores interiores, <math display="inline">\tilde{\epsilon }^{int}</math> , y, en la segunda etapa, se estiman los errores sobre aristas, <math display="inline">\tilde{\epsilon }^{edg}</math> . La estimación del error dual, denotada por <math display="inline">\tilde{\epsilon }</math> , se define como la suma de estas dos contribuciones, es decir, <math display="inline">\epsilon \approx \tilde{\epsilon }=\tilde{\epsilon }^{int}+\tilde{\epsilon }^{edg}</math> . Con esta estimación y la representación dada en la ecuación (9), se obtiene otra estimación del error en cantidades de interés, definida por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>E_4:=R^P\left(\tilde{\epsilon }\right)=R^P(\tilde{\epsilon }^{int})+R^P(\tilde{\epsilon }^{edg})\mbox{,}</math>
|}
| style="text-align: right;" | ( 17)
|}

donde <math display="inline">R^P\left(\tilde{\epsilon }^{int}\right)=\sum_{k=1}^{n_{elem}}\frac{R^D\left(\psi _k\right)R^P\left(\psi _k\right)}{a\left(\psi _k,\psi _k\right)}</math>  y <math display="inline">R^P\left(\tilde{\epsilon }^{edg}\right)=\sum_{l=1}^{n_{edg}}\frac{R^D\left(\chi _l\right)-a\left(\chi _l,\tilde{\epsilon }^{int}\right)}{a\left(\chi _l,\chi _l\right)}R^P\left(\chi _l\right)</math> .      

====6.1.1. Función burbuja====

En este trabajo, se utilizan diferentes funciones de tipo burbuja para determinar cuál de ellas ofrece la mejor estimación del error, tanto en energía como en cantidades de interés. Estas funciones se definen en un sistema de coordenadas <math display="inline">(\xi ,\eta )</math>  local. Luego se utilizan transformaciones de las usuales en el método del elemento finito para enviar dichas funciones a cada uno de los elementos de la malla o elementos globales. En esta ocasión, sobre las aristas, se parte de estas mismas funciones, pero transformadas de manera adecuada con rotaciones, traslaciones y dilataciones en el sistema <math display="inline">\left(\xi ,\eta \right)</math> , para que, al ser transformadas al elemento global, sean de soporte compacto sobre un subconjunto Λ<sub>l</sub>  que contenga la arista como el que se muestra en la [[#fig0020|figura 4]] . La burbuja sobre arista definida en el sistema local es enviada al elemento global Λ<sub>l</sub> .      

====6.1.2. Función burbuja escogida a posteriori (a partir de la solución de los problemas locales)====

Una de las opciones a la hora de escoger funciones de tipo burbuja es resolver problemas locales sobre cada elemento de la malla de elementos finitos, que es la estimación implícita descrita en el apartado 5. Con este proceso, se obtienen buenos resultados, aunque es costoso porque hay que resolver un problema por cada elemento y uno por cada arista.

====6.1.3. Función burbuja escogida a priori====

El primer tipo de función burbuja interior escogida a priori es como la presentada en [[#bib0180|[16]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\psi _x^1\left(\xi .\eta \right)=\psi _y^1\left(\xi ,\eta \right)=\left(1-\xi ^2\right)\left(1-\eta ^2\right)\mbox{,}</math>
|}
| style="text-align: right;" | 
|}

que tiene una derivada no nula en la dirección normal del contorno y un segundo tipo de función burbuja escogida a priori, de manera tal que tenga una derivada nula en el contorno:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\psi _x^2\left(\xi .\eta \right)=\psi _y^2\left(\xi ,\eta \right)=0,25(cos\left(\pi \xi \right)+1)(cos\left(\pi \eta \right)+1)\mbox{.}</math>
|}
| style="text-align: right;" | 
|}

==7. Resultados numéricos==

En este ejemplo, se tiene un dominio de material elástico lineal cuadrado, de dimensiones 1 × 1 m<sup>2</sup> , con un orificio cuadrado en el medio de 0,5 × 0,5 m<sup>2</sup> . El módulo de Young es E = 100 MPa; el coeficiente de Poisson es ν = 0,1 y la densidad es ρ = 1.000 kg/m<sup>3</sup> . Se imponen condiciones de Dirichlet homogéneas en el extremo izquierdo. Para el problema primal, se aplica una fuerza de 250 N, repartida en el segmento derecho de la parte superior del dominio, que resulta de aplicar una presión de 1.000 N/m sobre los 0,25 m del lado superior situados sobre la columna vertical. Las condiciones de contorno, la fuerza aplicada sobre el dominio y la deformada del problema primal (para la solución en los desplazamientos) se pueden observar en la [[#fig0025|figura 5]] . El problema se resolvió con el método de los elementos finitos, utilizando elementos cuadriláteros con funciones de forma bilineales, mediante Cast3 M [[#bib0200|[20]]] . La cantidad de interés escogida para el análisis es el desplazamiento del punto <math display="inline">\mbox{x}_0=(0,75;0,75)</math> . Para ello, se resuelve un problema auxiliar, que consiste en aplicar sobre el punto <math display="inline">\mbox{x}_0</math>  del dominio una fuerza <math display="inline">\mbox{j}</math>  de módulo 1 N. Utilizando el funcional dado en la ecuación (7), se escoge <math display="inline">\mbox{q}=\left(\frac{1}{\sqrt{2}},\frac{1}{\surd 2}\right)</math> . Las condiciones de contorno son iguales a las del problema primal. Las condiciones de contorno, la fuerza aplicada y la deformada (para la solución en los desplazamientos) del problema dual se pueden observar en la [[#fig0030|figura 6]] .

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_530383077-1-s2.0-S0213131515000504-gr5.jpg|center|284px|Deformada primal.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 5.

Deformada primal.

</span>
|}

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_530383077-1-s2.0-S0213131515000504-gr6.jpg|center|288px|Deformada dual.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 6.

Deformada dual.

</span>
|}

El problema se resolvió sobre una malla muy fina de 202.800 elementos, sobre la cual se obtuvo un valor de referencia para la norma energética de los desplazamientos de <math display="inline">\left\|\mbox{u}\right\|_E=3,7827</math>  y, para la cantidad de interés evaluada en dicha solución, de <math display="inline">J\left(\mbox{u}\right)=-5,45171\times 10^{-2}</math> . El problema se resolvió utilizando cuatro mallas más gruesas y, con el fin de no utilizar el espacio ocupado por la solución de referencia indicada anteriormente, en cada prueba se recalculó una nueva solución de referencia sobre la malla utilizada para la estimación. Los resultados pueden observarse en la [[#tbl0005|tabla 1]] , donde el número de elementos de la malla de cálculo y la utilizada para la estimación son denotados por n<sub>elem</sub>  y n<sub>eles</sub> , respectivamente; el error en la cantidad de interés, por J('''e''' ), y las representaciones, por <math display="inline">Rep_1=R^D\left(e\right)+R^P\left(\mbox{u}_H^\ast \right),Rep_2=R^P\left(\epsilon \right)+R^P\left(\mbox{u}_H^\ast \right)</math> , <math display="inline">R^P(\epsilon )</math>  y <math display="inline">R^D\left(e\right)</math> . El residuo primal evaluado en la solución del problema dual es denotado por <math display="inline">R^P\left(\mbox{u}_\mbox{H}^\ast \right)</math> ; la cantidad de interés evaluada en la solución, por <math display="inline">J\left(\mbox{u}_\mbox{H}\right)</math> ; la norma energética del error primal, de la solución, de la estimación implícita, de la estimación explícita con la burbuja <math display="inline">\psi ^1</math>  y de la estimación explícita con la burbuja <math display="inline">\psi ^2</math> , por <math display="inline">\left\|e\right\|_E</math> , <math display="inline">\left\|\mbox{u}_\mbox{H}\right\|_E</math> , <math display="inline">\left\|\overline{e}\right\|_E</math> , <math display="inline">\left\|e_1\right\|_E</math>  y <math display="inline">\left\|e_2\right\|_E</math> , respectivamente. Haciendo el análisis en energía en este problema, se observa que, con la estimación implícita y con la burbuja <math display="inline">\psi ^1</math> , se obtiene un índice de efectividad con una media del 81%, mientras que, con la burbuja <math display="inline">\psi ^2</math> , la media es del 61%. También se puede ver que, con una malla de 1.728 elementos, se obtiene un error relativo del 12%; es necesaria una malla de más de 120.000 elementos para disminuir este error a menos del 2%. Sin embargo, el propósito de este trabajo se centra en la cantidad de interés y, en este caso, el error relativo obtenido es del 2% con 1.728 elementos.

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 1.

Resultados de referencia sobre diferentes mallas y norma energética de los errores estimados

|-

! n<sub>elem</sub>
! n<sub>eles</sub>
! <math display="inline">J\left(e\right)=Rep_1=Rep_2</math>
! <math display="inline">J\left(\mbox{u}_\mbox{H}\right)</math>
! <math display="inline">\left|\left|e\right|\right|_E</math>
! <math display="inline">\left|\left|\mbox{u}_\mbox{H}\right|\right|_E</math>
! <math display="inline">R^D\left(e\right)=R^P\left(\epsilon \right)</math>
! <math display="inline">\mbox{R}^P\left(\mbox{u}_H^\ast \right)</math>
! <math display="inline">\left|\left|\overline{e}\right|\right|_E</math>
! <math display="inline">\left|\left|e_1\right|\right|_E</math>
! <math display="inline">\left|\left|e_2\right|\right|_E</math>
|-

! 108
| 1296
| -5,08071E-3
| -4,67E-2
| 1,1089
| 3,55
| -5,49E-3
| 4,17E-4
| 0,8854
| 0,8765
| 0,6698
|-

! 192
| 2304
| -3,66567E-3
| -4,88E-2
| 0,9279
| 3,62
| -3,95E-3
| 2,93E-4
| 0,7545
| 0,7466
| 0,5678
|-

! 432
| 5184
| -2,31288E-3
| -5,09E-2
| 0,7177
| 3,69
| -2,48E-3
| 1,75E-4
| 0,5878
| 0,5873
| 0,4341
|-

! 1728
| 20736
| -1,09475E-3
| -5,28E-2
| 0,4652
| 3,74
| -1,16E-3
| 0,74E-4
| 0,3770
| 0,3877
| 0,290
|}

En la estimación del error, se ensayó con tres tipos diferentes de funciones de tipo burbuja para observar la influencia de estas sobre los resultados.

'''Ensayo 1.'''  Sobre cada uno de los elementos de una malla gruesa, se resuelve un problema con condiciones de Dirichlet homogéneas, con el residuo como término fuente, y se escoge como burbuja la solución de este problema (una burbuja para cada elemento). Las estimaciones del error en la cantidad de interés dadas en las ecuaciones (14) y (15), <math display="inline">\mbox{E}_{\alpha }</math> , el error relativo <math display="inline">\eta _{\alpha }=\left|\frac{E_{\alpha }}{J\left(\mbox{u}\right)}\right|</math>  y los índices de efectividad <math display="inline">\Theta _{\alpha }=\frac{E_{\alpha }}{J(e)}\ast 100\mbox{ }\%,\mbox{ }\mbox{ para}\mbox{ }\alpha =1,2</math> , se pueden observar en la [[#tbl0010|tabla 2]] . La media de los índices de efectividad en este ensayo es del 98,4% y se obtiene convergencia del orden de 1,2. Sobre la malla de 1.728 elementos, el error relativo estimado es inferior al 2%. En este ensayo, se resuelven n<sub>elem</sub>  problemas sobre mallas de 12 elementos cada una y n<sub>edg</sub>  problemas sobre mallas de 12 elementos cada una, para un total de (n<sub>elem</sub> + n<sub>edg</sub> ) problemas. Para la malla de 2.700 elementos y 2.880 nodos, esto representa un poco más de 5.580 problemas resueltos.

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 2.

Estimación implícita sobre mallas uniformes

|-

! n<sub>elem</sub>
! E<sub>1</sub>
! E<sub>2</sub>
! <math display="inline">\eta _1</math>
! <math display="inline">\eta _2</math>
! <math display="inline">\Theta _1</math>
! <math display="inline">\Theta _2</math>
|-

! 108
| -5,00915E-3
| -5,04821E-3
| 9,66E-02
| 9,74E-02
| 98,59
| 99,36
|-

! 192
| -3,70759E-3
| -3,7284E-3
| 7,05E-02
| 7,09E-02
| 101,1
| 101,7
|-

! 432
| -2,31759E-3
| -2,32338E-3
| 4,35E-02
| 4,36E-02
| 100,2
| 100,4
|-

! 1728
| -1,01868E-3
| -1,01488E-3
| 1,88E-02
| 1,88E-02
| 93,05
| 92,70
|}

En los dos ensayos siguientes, se resolvió el problema sobre cuatro mallas uniformes, con submallas no estructuradas para calcular las estimaciones <math display="inline">\tilde{e}_j</math>  y <math display="inline">\tilde{\epsilon }_j</math>  del error primal y el dual, con las burbujas <math display="inline">\psi ^j</math> , para j = 1, 2. Los resultados se pueden observar en las [[#tbl0015|Tabla 3]]  and [[#tbl0020|Tabla 4]] , donde las estimaciones son <math display="inline">E_{1j}=R^D\left(\tilde{e}_j\right)</math>  y <math display="inline">E_{2j}=R^P\left(\tilde{\epsilon }_j\right)</math> , los errores relativos son definidos por <math display="inline">\eta _{1j}=\left|\frac{E_{1j}}{J(\mbox{u})}\right|</math>  y <math display="inline">\eta _{2j}=\left|\frac{E_{2j}}{J(\mbox{u})}\right|</math> , y los índices de efectividad son definidos por <math display="inline">\Theta _{1j}=\frac{E_{1j}}{J(e)}\ast 100\%</math>  y <math display="inline">\Theta _{2j}=\frac{E_{2j}}{J(e)}\ast 100\%</math> , para j = 1,2.

<span id='tbl0015'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 3.

Estimación sobre mallas uniformes con <math display="inline">\psi ^1\left(\xi ,\eta \right)</math>

|-

! n<sub>elem</sub>
! <math display="inline">\mbox{E}_{11}</math>
! <math display="inline">\mbox{E}_{21}</math>
! <math display="inline">\eta _{11}</math>
! <math display="inline">\eta _{21}</math>
! <math display="inline">\Theta _{11}</math>
! <math display="inline">\Theta _{21}</math>
|-

! 108
| -3,10714E-3
| -3,20829E-3
| 5,99E-2
| 6,19E-02
| 61,1
| 63,1
|-

! 192
| -2,26008E-3
| -2,34881E-3
| 4,30E-02
| 4,46E-02
| 61,6
| 64,0
|-

! 432
| -1,39022E-3
| -1,45637E-3
| 2,61E-02
| 2,73E-02
| 60,1
| 62,9
|-

! 1728
| -5,97947E-4
| -6,4447E-4
| 1,10E-02
| 1,19E-02
| 54,6
| 58,8
|}

<span id='tbl0020'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 4.

Estimación sobre mallas uniformes con <math display="inline">\psi ^2\left(\xi ,\eta \right)</math>

|-

! n<sub>elem</sub>
! <math display="inline">\mbox{E}_{12}</math>
! <math display="inline">\mbox{E}_{22}</math>
! <math display="inline">\eta _{12}</math>
! <math display="inline">\eta _{22}</math>
! <math display="inline">\Theta _{12}</math>
! <math display="inline">\Theta _{22}</math>
|-

! 108
| -2,42236E-3
| -2,45699E-3
| 4,67E-02
| 4,74E-02
| 47,6
| 48,3
|-

! 192
| -1,75857E-3
| -1,78471E-3
| 3,34E-02
| 3,39E-02
| 47,9
| 48,6
|-

! 432
| -1,08154E-3
| -1,10251E-3
| 2,03E-02
| 2,07E-02
| 46,7
| 47,6
|-

! 1.728
| -4,6204E-4
| -4,76812E-4
| 8,56E-03
| 8,83E-03
| 43,5
| 43,5
|}

'''Ensayo 2.'''  Se escogió la burbuja interior, definida por <math display="inline">\psi ^1=\lbrace \psi _x^1,\left(1+\nu \right)\psi _y^1\rbrace </math> , donde <math display="inline">\psi _x^1=\psi _y^1=\left(1-\xi ^2\right)\left(1-\eta ^2\right)</math> . Con esta elección, se obtuvieron los resultados que se pueden observar en la [[#tbl0015|tabla 3]] .      

En este ensayo, la media de los índices de efectividad para las estimaciones en cantidad de interés es del 60,8% y se obtiene convergencia de orden 1. Sobre la malla de 1.728 elementos, el error relativo de las estimaciones es inferior al 2%.

'''Ensayo 3.'''  Se escogió la burbuja interior definida por <math display="inline">\psi ^2\left(\xi ,\eta \right)=\lbrace \psi _x^2,\left(1+\nu \right)\psi _y^2\rbrace </math> , donde <math display="inline">\psi _x^2=\psi _y^2=\frac{1}{4}\left(cos\left(\pi \xi \right)+1\right)\left(cos\left(\pi \eta \right)+1\right)</math> . Con esta elección, se obtuvieron los resultados que se pueden observar en la [[#tbl0020|tabla 4]] .      

La media de los índices de efectividad en este ensayo es del 46,7% y se obtiene convergencia del orden de 1,2. Sobre la malla de 1.728 elementos, el error relativo de las estimaciones es inferior al 1%.

En la [[#fig0035|figura 7]] , se puede observar el mapa de los errores estimados con <math display="inline">E_2</math> . Para dibujarlo, se procedió de la forma siguiente. Los errores interiores se asignaron por elementos tal como se construyeron. Los errores sobre aristas se distribuyeron de manera que, si el subdominio utilizado se interseca con dos elementos, la contribución se reparte a partes iguales para cada uno de los elementos involucrados, mientras que, si el subdominio está contenido en un solo elemento, como es el caso de las aristas frontera, la contribución se asigna completamente a dicho elemento. Los errores son mayores en las esquinas interiores, donde surgen singularidades, producto de las tensiones producidas por las fuerzas aplicadas y la geometría del dominio. Para analizar estos sectores, se podría utilizar un método de refinamiento adaptativo, utilizando las estimaciones presentadas aquí, para colocar un mayor número de elementos en estas zonas y tratar de disminuir el error de la solución aproximada de la cantidad de interés; sin embargo, no se ha procedido así en este trabajo, porque estamos pensando en utilizar estas estimaciones en problemas transitorios. En la [[#fig0040|figura 8]]  se puede ver en escala logarítmica, el error relativo de la solución, cuando el problema es resuelto utilizando mallas uniformes, donde <math display="inline">J(e)</math>  representa el error en la cantidad de interés. El orden de convergencia de cada método tiene la misma magnitud que el error en la cantidad de interés de referencia, que en este problema es de 1.2.

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_530383077-1-s2.0-S0213131515000504-gr7.jpg|center|328px|Mapa de errores estimados.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 7.

Mapa de errores estimados.

</span>
|}

<span id='fig0040'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_530383077-1-s2.0-S0213131515000504-gr8.jpg|center|348px|Error relativo de las estimaciones.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 8.

Error relativo de las estimaciones.

</span>
|}

El índice de efectividad de las estimaciones analizadas se puede observar en la [[#fig0045|figura 9]] , donde la notación es la indicada al comienzo del apartado. En dicha figura, puede apreciarse que la norma en energía del error estimado implícitamente tiene un índice de efectividad muy similar al obtenido con la estimación con la burbuja <math display="inline">\psi ^1</math> . También puede observarse que, para la cantidad de interés, es la estimación implícita la que presenta mejores índices de efectividad, seguida de la estimación con la burbuja <math display="inline">\psi ^1</math>  y, por último, de la estimación con la burbuja <math display="inline">\psi ^2</math> . Aunque los mejores resultados se obtienen con la estimación implícita, esta tiene el costo computacional más elevado, pues se resuelven <math display="inline">n_{elem}+n_{edg}</math>  problemas de contorno mientras que, al estimar el error utilizando funciones burbuja, solo se realizan productos de matrices con vectores sobre cada elemento.

<span id='fig0045'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_530383077-1-s2.0-S0213131515000504-gr9.jpg|center|345px|Índice de efectividad de las estimaciones.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 9.

Índice de efectividad de las estimaciones.

</span>
|}

==8. Conclusiones==

Se ha estimado el error ensayando con tres funciones de tipo burbuja, definidas con soporte compacto. La primera es implícita, escogida como la combinación ortogonal de dos funciones, en que una es la solución de un problema de Dirichlet sobre cada elemento de la malla y la otra es la solución de un problema sobre un subdominio que contiene cada arista de la malla. La segunda y la tercera fueron escogidas a priori, y una de las diferencias entre ellas es que las derivadas normales en el contorno del dominio son nulas para una y no nulas para la otra. Entre los tres tipos de funciones, las implícitas dieron mejores índices de efectividad, aunque el orden de convergencia fue similar con todas las estimaciones.

El problema analizado presenta singularidades debidas a las esquinas reentrantes presentes en el dominio, lo cual hace necesario un número grande de elementos para obtener una buena aproximación. Aunque, en energía, este hecho es muy evidente, pues se necesitan cerca de 20.000 elementos para alcanzar el régimen de convergencia asintótico, en el análisis de la cantidad de interés se requieren menos elementos.

A medida que el número de elementos de las mallas se incrementa, los errores relativos de las estimaciones decrecen, lo cual es deseable en un buen estimador. En el error estimado en la cantidad de interés, utilizando implícitamente y explícitamente diferentes tipos de funciones burbuja en el problema mecánico analizado, se obtuvo convergencia del orden de 1,2.

Para la cantidad de interés, el índice de efectividad medido con la estimación implícita tiene un media del 98,4%; con la burbuja <math display="inline">\psi ^1</math> , del 60,8% y, con la burbuja <math display="inline">\psi ^2</math> , del 46,7%.      

Se ha analizado el efecto de proyectar una solución desde una malla gruesa de cuadriláteros y funciones de forma bilineales sobre una malla no estructurada y se ha observado que, en esta situación, aparece un nuevo término en las diferentes representaciones del error en cantidades de interés, el cual es fácil de calcular, una vez obtenida una solución del problema dual. Este término afecta el índice de efectividad, pero no el orden de convergencia, y es nulo cuando las mallas pertenecen al mismo espacio.

Una vez obtenidos los resultados del problema analizado, se puede sugerir el uso del método semiexplícito para estimar el error en cantidades de interés en los problemas mecánicos, en aquellos casos en que se requiera una buena precisión, o el método explícito utilizando funciones de tipo burbuja, si se quiere una mayor velocidad de cálculo, aunque sea a costa de perder un poco de precisión. Con el método explícito, no se requiere resolver problemas locales sobre todo el dominio, lo cual permite reducir el tiempo de cálculo, pues solo hay que realizar multiplicaciones de matrices y vectores. Puede ser de gran utilidad para los problemas de adaptatividad o para los problemas transitorios, donde es preciso estimar el error en múltiples ocasiones.

==Agradecimientos==

Este trabajo ha contado con el apoyo del Consejo de Desarrollo Científico, Humanístico, Tecnológico y de las Artes (CDCHTA) de la Universidad de Los Andes, a través del proyecto I -1204-09-05-C.

==Referencias==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0105'></span>
[[#bib0105|[1]]] M. Ainsworth, J. Oden; A posteriori error estimation in finite element analysis; Wiley Inter-Science (2000)</li>
<li><span id='bib0110'></span>
[[#bib0110|[2]]] P. Díez, N. Parés, A. Huerta; Recovering lower bounds of the error by postprocessing implicit residual a error estimates; International Journal for Numerical Methods in Engineering, 56 (10) (2003), pp. 1465–1488</li>
<li><span id='bib0115'></span>
[[#bib0115|[3]]] R. Cottereau, P. Díez, A. Huerta; Strict Error Bounds for Linear Solid Mechanics Problems using a Subdomain-Based Flux-Free Method; Computational Mechanics, 44 (4) (2009), pp. 533–547</li>
<li><span id='bib0120'></span>
[[#bib0120|[4]]] N. Parés, P. Díez, A. Huerta; Exact bounds for linear outputs of the advection-diffusion-reaction equation using flux-free error estimates; SIAM Journal of Scientific Computing, 31 (4) (2009), pp. 3064–3089</li>
<li><span id='bib0125'></span>
[[#bib0125|[5]]] D. Aubry, D. Lucas, B. Tie; Adaptive strategy for transient/coupled problems applications to thermoelasticity and elastodynamics; Comput. Methods Appl. Mech. Engrg., 176 (1999), pp. 41–50</li>
<li><span id='bib0130'></span>
[[#bib0130|[6]]] W. Nils-Erik, X. Li; Adaptive finite element procedures for linear and non-linear dynamics; International Journal for numerical methods in engineering, 46 (1999), pp. 1781–1802</li>
<li><span id='bib0135'></span>
[[#bib0135|[7]]] G. Calderón, P. Díez; Análisis de diferentes estimadores de error de postproceso para adaptatividad orientada al resultado; Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería, 22 (2006), pp. 193–214</li>
<li><span id='bib0140'></span>
[[#bib0140|[8]]] W. Bangerth, R. Rannacher; Adaptive finite element methods for differential equations; Birkhäuser Verlag (2003)</li>
<li><span id='bib0145'></span>
[[#bib0145|[9]]] M. Botkin, H. Wang; An adaptive mesh refinement of quadrilateral finite element meshes based upon a posteriori error estimation of quantities of interest: linear static response; Engineering with Computers, 20 (2004), pp. 31–37</li>
<li><span id='bib0150'></span>
[[#bib0150|[10]]]  Zienkiewicz,  Zhu; The superconvergent patch recovery (SPR) and adaptive finite element refinement; Computer Methods in Applied Mechanics and Engineering., 101 (1992), pp. 207–224</li>
<li><span id='bib0155'></span>
[[#bib0155|[11]]] J.J. Ródenas, O.A. González-Estrada, P. Díez, F.J. Fuenmayor; Accurate recovery-based upper error bounds for the extended finite element framework; Computer Methods in Applied Mechanics and Engineering, 199 (37–40) (2010), pp. 2607–2621</li>
<li><span id='bib0160'></span>
[[#bib0160|[12]]] P. Díez, J.J. Ródenas, O.C. Zienkiewicz; Equilibrated Patch Recovery error estimates: simple and accurate upper bounds of the error; International Journal for Numerical Methods in Engineering, 69 (2007), pp. 2075–2098</li>
<li><span id='bib0165'></span>
[[#bib0165|[13]]] F. Verdugo, P. Díez; Computable bounds of functional outputs in linear visco-elastodynamics; Computer Methods in Applied Mechanics and Engineering, 245 (2012), pp. 313–330</li>
<li><span id='bib0170'></span>
[[#bib0170|[14]]] P. Díez, J.J. Egozcue, A. Huerta; A posteriori error estimation for standard finite element analysis; Computer Methods in Applied Mechanics and Engineering., 163 (1998), pp. 141–157</li>
<li><span id='bib0175'></span>
[[#bib0175|[15]]] E. Stein, M. Rüter, S. Ohnimus; Error-controlled adaptive goal-oriented modeling and finite element approximations in elasticity; Comput. Methods Appl. Mech. Engrg., 196 (2007), pp. 3598–3613</li>
<li><span id='bib0180'></span>
[[#bib0180|[16]]] R. Rosales, P. Díez; Estima de error residual explícita para cantidades de interés utilizando funciones burbuja; Revista Internacional Métodos Numéricos para Cálculo y Diseño en Ingeniería., 25 (2009), pp. 337–357</li>
<li><span id='bib0185'></span>
[[#bib0185|[17]]] S. Prudhomme, J.T. Oden; On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors; Comput. Methods Appl. Mech., 176 (1999), pp. 313–331</li>
<li><span id='bib0190'></span>
[[#bib0190|[18]]] W. Bangerth, R. Rannacher; Adaptive finite element methods for differential equations; Birkhäuser Verlag, Basel–Boston–Berlin (2003)</li>
<li><span id='bib0195'></span>
[[#bib0195|[19]]] J.T. Oden, S. Prudhomme; Goal-oriented error estimation and adaptivity for the finite element method; Computers and Mathematics with Applications, 41 (2001), pp. 735–756</li>
<li><span id='bib0200'></span>
[[#bib0200|[20]]] http://www-cast3m.cea.fr</li>
</ol>