Upload
trinhkien
View
212
Download
0
Embed Size (px)
Citation preview
Thiago Rafael Bremm
Simulação física de corpos deformáveis com técnicasem malha: implementação de método baseado em
pontos.
Florianópolis, Santa Catarina
Novembro de 2012
Thiago Rafael Bremm
Simulação física de corpos deformáveis com técnicasem malha: implementação de método baseado em
pontos.
Trabalho de conclusão de curso apresentadocomo parte dos requisitos para obtenção do graude Bacharel em Ciências da Computação
Orientador:
Prof. Dr. rer.nat. Aldo von Wangenheim
Co-Orientadores Tiago de Holanda Cunha Nobrega, Jeferson Ramos
Florianópolis, Santa Catarina
Novembro de 2012
Resumo
Trabalho que realiza pesquisa de modelos físicos e matemáticos que descrevem a dinâmicados sólidos e os métodos computacionais para sua simulação, expondo a divisão em duasgrandes classes. Estudando os clássicos com uso de malha, bem estabelecidos e amplamenteaplicados desde os anos 60, em comparação aos sem malha de pesquisa recente. Explorandoo potencial desta nova abordagem através da implementação da técnica de animação elásticaplástica baseada em pontos, verificando seus benefícios e dificuldades com a construção de umsimulador interativo para aplicações de tempo real.
Contents
List of Figures
1 Introdução p. 8
1.1 Objetivos Gerais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 9
1.2 Objetivos Específicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 9
1.3 Requisitos de sistema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 10
1.3.1 Funcionais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 10
1.3.2 Não funcionais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 10
2 Revisão Sistemática da Literatura p. 12
2.1 Protocolo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 12
2.2 Parâmetros de exclusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 13
2.3 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 14
2.3.1 Mecânica Contínua, dinâmica dos sólidos, elasticidade e plasticidade p. 14
2.3.2 Simulação de corpos deformáveis em tempo real . . . . . . . . . . . p. 15
2.3.3 Métodos com malha e sem malha . . . . . . . . . . . . . . . . . . . p. 15
2.3.4 Dos métodos com malha . . . . . . . . . . . . . . . . . . . . . . . . p. 15
2.3.5 Dos métodos livres de malha . . . . . . . . . . . . . . . . . . . . . . p. 16
3 Mecânica Contínua p. 17
3.1 Deslocamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 17
3.2 Gradiente do Deslocamento . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 18
3.3 Tensão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 19
3.4 Estresse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 20
4 Modelos Computacionais p. 23
4.1 Técnicas Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 23
4.1.1 Massa-Mola . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 24
4.1.2 FEM - Método dos Elementos Finitos . . . . . . . . . . . . . . . . . p. 28
4.1.3 Construção e preservação da Malha . . . . . . . . . . . . . . . . . . p. 30
4.2 Técnicas Meshless . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 31
5 Discretização do contínuo p. 35
5.1 Método de Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 35
5.2 Amostragem e função aproximadora . . . . . . . . . . . . . . . . . . . . . . p. 36
5.3 Least Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 36
5.4 Weighted Least Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 39
5.4.1 Função de peso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 39
5.4.2 Partição do domínio . . . . . . . . . . . . . . . . . . . . . . . . . . p. 41
5.5 Moving Least Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 42
6 Simulação elástica e plástica baseada em pontos p. 44
6.1 Visão geral do método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 44
6.2 Inicialização dos phyxels . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 46
6.3 Calculando o gradiente do deslocamento . . . . . . . . . . . . . . . . . . . . p. 46
6.3.1 Aplicando MLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 46
6.3.2 Inversão por SVD . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 48
6.4 Calculando a tensão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 50
6.4.1 Tensão elástica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 50
6.4.2 Tensão plástica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 51
6.5 Calculando o estresse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 52
6.6 Calculando as forças resultantes . . . . . . . . . . . . . . . . . . . . . . . . p. 53
6.7 Integração temporal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 54
6.8 Visualização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 55
6.8.1 Atualização dos surfels pelo deslocamento . . . . . . . . . . . . . . p. 56
7 Implementação p. 58
7.1 Recursos e bibliotecas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 58
7.2 Arquitetura do protótipo . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 59
7.3 Carregamento de objetos . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 59
7.4 Simulação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 61
7.4.1 Configuração dos Phyxels . . . . . . . . . . . . . . . . . . . . . . . p. 62
7.4.2 Aplicando SVD no cálculo da inversa de A . . . . . . . . . . . . . . p. 62
7.4.3 Colisão com superfícies estáticas . . . . . . . . . . . . . . . . . . . . p. 63
7.4.4 Aplicação de forças externas . . . . . . . . . . . . . . . . . . . . . . p. 64
7.4.5 Atualização da visualização . . . . . . . . . . . . . . . . . . . . . . p. 64
7.5 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 65
8 Conclusão p. 67
Referências Bibliográficas p. 69
List of Figures
3.1 Variação da distância entre dois pontos de um contínuo após sofrer deformação. p. 19
4.1 Malha definida por polígonos com seus vértices e arestas . . . . . . . . . . . p. 23
4.2 Malhas de polígonos representando corpos tridimensionais. (CLKER, ) . . . p. 24
4.3 Compressão sofrida por uma mola em decorrência da aplicação de uma força
fi,com um deslocamento equivalente a L−L′. . . . . . . . . . . . . . . . . . p. 25
4.4 Exemplo ilustrando massas conectadas entre si por molas. . . . . . . . . . . p. 26
4.5 Simulação em tempo real de tecido usando massa-mola. . . . . . . . . . . . . p. 27
4.6 Elemento 2D sofrendo deformação e as forças elásticas atuando nos nodos. . p. 28
4.7 Malha em sistema de massa-mola planejada para comportamento isotrópico
(NEALEN et al., 2006) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 31
4.8 Malha representando estruturas anatômicas de um tronco humano, com a dis-
tribuição dos ângulos dos poliedros. (FOTEINOS; CHRISOCHOIDES, 2011) p. 32
4.9 a) Malha com buracos. b) Fidelidade da superfície aprimorada com técnica
de rejeição de pontos. c) Corte transversal mostrando a estrutura interna da
malha. (FOTEINOS; CHRISOCHOIDES, 2011) . . . . . . . . . . . . . . . p. 32
4.10 Simulação de rachadura com EFGM, onde o domínio é representado apenas
por pontos, sem necessidade de ligações. (BELYTSCHKO; TABBARA, 1996) p. 33
4.11 Dinâmica das partículas em protótipo de simulação de fluídos baseado em
SPH. (NOBREGA, 2010) . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 34
5.1 Dois exemplos de superfícies geradas a partir de duas amostras de pontos
aproximadas pelos polinômios gerados por mínimos quadrados. (NEALEN;
DARMSTADT, 2004) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 38
5.2 A curva em azul representa uma aproximação por mínimos quadrados, onde
a influência dos pontos é equivalente e global, o polinômio gerado tende a
se aproximar de todos com o menor erro possível. A curva rosada representa
uma aproximadora gerada por mínimos quadrados ponderados, onde a distân-
cia dos pontos ao centro de análise é considerada. Desta forma a influência
de pontos mais afastados sobre a curva é inferior à dos pontos próximos a x. . p. 40
5.3 Em gráfico o comportamento de três diferentes funções de peso com raios de
tamanho 1. A gaussiana em roxo (NEALEN; DARMSTADT, 2004), poly6
em verde (MÜLLER et al., 2004) e spiky kernel em amarelo (NOBREGA,
2010). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 41
5.4 Demonstração de partição de domínio com pontos de análise com abrangên-
cia de raios grandes, buscando cobrir toda a área do domínio. . . . . . . . . . p. 42
6.1 Demonstração do resultado da simulação de plasticidade e elasticidade com
a técnica baseada em pontos, que podem ser vistos na figura da esquerda.
(MÜLLER et al., 2004) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 45
6.2 Distribuição coplanar de partículas. . . . . . . . . . . . . . . . . . . . . . . . p. 50
7.1 Visão superficial dos componentes que compõe o C3DE. (SILVA et al., 2009) p. 59
7.2 Demonstração de resultado do algoritmo de preenchimento com phyxels. A
primeira linha, em destaque, exibe a contagem de colisões desta passada
através da figura, com os phyxels colocados nos pontos ondem o número
é impar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 61
7.3 Do modelo em malha de triângulos visto na esquerda, obteu-se a nuvem de
partículas exibida na direita em resultado do algoritmo de preenchimento. . . p. 61
7.4 Resultado da simulação de colisão entre modelo deformável e uma esfera
estática. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 63
7.5 Três quadros que exibem o resultado da ferramenta de pinçar, a caixa ver-
melha marca o local do seletor. . . . . . . . . . . . . . . . . . . . . . . . . p. 64
8
1 Introdução
Simular características do mundo real em ambiente computacional interessa pesquisadores
de diversas áreas. Suas aplicações exploram motivações práticas e teóricas que variam desde
entretenimento a modelos de predição, testes virtuais, prototipação de mecanismos, resistência
de materiais, reações químicas e fenômenos físicos em geral. Reproduzir tais efeitos com el-
evado grau de fidelidade representa um desafio para a computação, pela alta complexidade da
maioria dos modelos e o grande volume de dados envolvidos.
A priori de qualquer representação, deve-se compor uma abstração do fenômeno estudado
de forma a gerar uma aproximação razoável sob a ótica de certos graus de aceitação, que se
encaixe nos limites das máquinas e técnicas contemporâneas. Com o exemplo da dinâmica de
corpos rígidos, que descarta certas propriedades físicas da matéria para exibir o comportamento
geral do movimento, modelando com maior eficiência os resultados pertinentes da interação
entre objetos com suposta rigidez infinita. O que é bastante adequado para muitos casos de
uso, como jogos eletrônicos com apelo baseado em realismo, sistemas de partículas, e diversas
aplicações de engenharia.
Descartada a deformação ou outra alteração física do material, o que resta para o modelo
considerar é bem conhecido e dominado pela mecânica clássica. São questões relacionadas a
posição, velocidade, aceleração, força, momento linear e angular, energia cinética, torque, etc.
Computacionalmente falando, os desafios compõe-se em dominar situações com grande número
de objetos, onde se deve integrar o movimento resultante de cada, detectar e tratar rapidamente
eventuais colisões entre estes (SOUZA, 2011). Esta é uma abstração afastada da realidade de
muitos problemas, onde não se pode aceitar que corpos sejam completamente rígidos, pois
sua deformação influencia gravemente nos resultados. Como a reação de estruturas a fortes
pressões, representação virtual de tecido vivo, ou qualquer material maleável.
Passa-se a abordar o estudo de mecânica dos meios contínuos, área da física onde pode-se
encontrar a formalização dos fenômenos desta natureza. Considerando elasticidade, plastici-
dade, fluidez, fratura, derretimento, termodinâmica, entre outros não contemplados por este
9
trabalho (BOWEN, 2010).
Aceitar a mudança estrutural nas colisões macroscópicas com rigidez finita requerer uma
abstração que considere mais do que o movimento de uma forma geométrica no espaço. Surge
na questão o conceito de configuração de um corpo, onde o mesmo precisa ter sua forma con-
cisamente descrita, e passível de manipulação a cada quadro do evento. Isto não pode ser feito
de maneira leviana, pois a quantidade de informação e complexidade dos fenômenos envolvidos
pode facilmente saturar o sistema ou pecar na produção de resultados realistas.
Considerando isso, vários modelos para simulação de corpos deformáveis foram e estão
sendo desenvolvidos por vários pesquisadores. Ressaltando-se duas classes: aqueles que estru-
turam o modelo com malha (mesh), e os livres de malha (meshless). Sendo a primeira bem con-
ceituada e amplamente empregada desde a primeira metade do século XX (FELIPPA, 2001),
em oposição à segunda, cujo desenvolvimento e pesquisa são relativamente recentes (IDEL-
SOHN; ONATE, 2006). Esta monografia elabora um comparativo entre ambas e segue para o
enfoque do estudo de um método livre de malha baseado em pontos, desenvolvido no trabalho
de Müller et al. (2004). Buscando atingir resultados com performance adequada para aplicações
interativas, e avaliar o potencial desta classe de métodos.
1.1 Objetivos Gerais
Através deste trabalho, espera-se aglomerar conhecimentos sobre simulação física de cor-
pos deformáveis. Especializando no estudo e desenvolvimento de técnicas livres de malha,
disponíveis na literatura. Avaliando seu potencial de aplicação por meio da implementação e
prototipação de uma destas, escolhida com preocupação em realismo, projetando a futuro uso
em aplicações de caráter científico e educacional.
1.2 Objetivos Específicos
• Desenvolver componentes de biblioteca para simulação elástica e plástica com método
baseado em pontos.
• Criar protótipo de simulação em tempo real que responda à interação do usuário.
• Validar resultados visuais comparando com comportamento real.
• Avaliar desempenho e estabilidade.
10
• Expor os pontos sensíveis da técnica, relatando suas vantagens e deficiências.
• Avaliar potencial de aplicação em simulações cirúrgicas.
1.3 Requisitos de sistema
A elaboração do protótipo focou o apoio à validação visual do método com controle e
configurabilidade da simulação para caráter de teste, ainda com avaliação de apelo visual e
viabilidade de interação em tempo real com o usuário.
1.3.1 Funcionais
• O simulador deve carregar arquivos no formato OBJ.
• O simulador deve permitir que o usuário configure a construção do modelo do objeto
carregado: granularidade e número de vizinhos.
• O simulador deve permitir que o usuário configure todos os parâmetros do material.
• O simulador deve permitir que os parâmetros sejam mudados ao longo da simulação.
• O simulador deve oferecer visão interativa, com movimento de câmera.
• O simulador deve oferecer ferramenta de interação com os objetos da cena com uso do
mouse: pinçar e tocar.
• A simulação deve ser controlável, permitindo reinício, pausa, continuidade, parada e
avanço por unidade de um intervalo de tempo.
• A simulação deve ser em tempo real.
• Ao longo da simulação, o usuário poderá escolher entre ver a superfície opaca do objeto,
ou transparente exibindo as partículas interiores.
1.3.2 Não funcionais
• O código será feito na linguagem C++.
• O programa e componentes de biblioteca devem ser compiláveis em ambiente Linux e
Windows.
11
• O código com suas classes e funções deve se encaixar nos padrões da biblioteca C3DE
do grupo INCoD.
12
2 Revisão Sistemática da Literatura
Este capítulo apresenta o desenvolvimento e a aplicação de uma revisão sistemática da lit-
eratura acerca de simulação física computadorizada de corpos deformáveis. Seguindo os passos
propostos por Kitchenham (2004) para a produção de uma revisão imparcial e reproduzível:
identificação da necessidade de uma revisão, desenvolvimento de um protocolo, identificação
da pesquisa, seleção de estudos de interesse, validação da qualidade, extração de conhecimento,
sintetização de dados e descrição da revisão.
2.1 Protocolo
A revisão sistemática foi dirigida pela formulação de um protocolo que segue o padrão
definido por Kitchenham (2004).
Uma série de perguntas iniciais definem o escopo dos termos chave para a pesquisa:
• Que fenômenos regem a deformação de corpos?
• Quais os métodos computacionais para simular esse efeito?
• Quais as áreas de aplicação?
• Quais técnicas são mais utilizadas em cada área de aplicação?
• Como se modela os objetos para a simulação?
• Como funcionam as técnicas com malha?
• Como funcionam as técnicas sem malha?
• Que vantagens uma possui sobre a outra?
• Qual a viabilidade de simular em tempo real com técnicas sem malha?
• Que etapas envolvem o processo de simular um corpo deformável?
13
• Como tratar colisão?
• Como realizar a integração do movimento?
Os seguintes conjuntos de termos chave para a busca foram definidos a partir das questões
levantadas:
• soft body dynamics, continuum mechanics, elasticity, plasticity
• real time softbody simulation
• mesh technics, meshless technics, mesh free
• tridimensional, object, sampling, representation, volume
• collision detection
• movement integration
A busca pelos conjuntos acima foi realizada nos seguintes motores de busca:
• Google Scholar (http://scholar.google.com)
• IEEExplore (http://ieeexplore.ieee.org/Xplore/guesthome.jsp)
• ACM Digital Library (http://portal.acm.org/)
• Science Direct (http://www.sciencedirect.com)
• Google (http://google.com)
2.2 Parâmetros de exclusão
Para refinar os resultados iniciais da pesquisa foram definidos critérios para exclusão de
documentos, cujo conteúdo é estimado a partir de título, resumo e posição nos motores de
busca.
Para verificar relevância, a análise sobre os seguintes itens foi feita:
• Título que não envolva o escopo deste trabalho.
• Resumo exibindo assunto incoerente com os aqui estudados.
14
• Títulos e resumos com tema repetitivo aos demais.
• Resumo pouco objetivo ou mal escrito.
• Como os motores de busca ordenam por maior afinidade com os termos usados, entre
outras características de relevância, tomou-se como limite máximo os trinta primeiros
itens retornados.
Após aplicação dos primeiros parâmetros de exclusão, a seleção continua com a leitura
preliminar do conteúdo dos documentos, em especial das informações contidas nas introduções
e capítulos iniciais, adaptando os critérios de exclusão sobre esses. Mais uma vez avaliando a
coesão entre os assuntos relatados e os de interesse, bem como clareza e qualidade expressadas
nos textos.
2.3 Resultados
As amostragens de artigos obtidos nas buscas efetuadas, com filtragem pelos parâmetros de
exclusão em análises superfíciais dos conteúdos, convergiram para 53 documentos. Os quais
foram separados em classes de acordo com as chaves de busca das quais resultaram, e uma
nova classificação foi feita avaliando seu conteúdo, especialmente introdução, estrutura de tópi-
cos e conclusão. Pela nova seleção, restaram 16 artefatos cujo conteúdo foi utilizado para a
elaboração deste trabalho. A seguir, a extração por classes apresentando breve resumo dos
documentos mais importantes.
2.3.1 Mecânica Contínua, dinâmica dos sólidos, elasticidade e plastici-dade
Devido ao caráter teórico fundamental, este fragmento da busca teve ainda outro critério,
priorizando a seleção de livros bem conceituados ao invés de artigos ou trabalhos específicos.
Dos livros encontrados, optou-se por aqueles de caráter introdutório permitisse compreen-
são avançada destas teorias. Neste quesito destaca-se o livro Introduction to Continuum Me-
chanics for Engineers de Bowen (2010), no qual pode-se encontrar as definições básicas e
elaboração das teorias e modelos matemáticos que governam a física dos corpos deformáveis,
incluindo revisões de cálculo e álgebra vetorial. O conteúdo deste livro permite a compreensão
e confere capacidade crítica para estudar os modelos computacionais que simulam tais fenô-
menos.
15
Em complemento ao assunto abordado pelo livro anteriormente descrito, a obra Introduc-
tion to Tensor Calculus and Contin-uum Mechanics de Heinbockel (2001) reforça a teoria da
dinâmica dos sólidos, e ainda dá ênfase especial ao assunto de tensores, ferramenta matemática
muito importante na mecânica contínua.
2.3.2 Simulação de corpos deformáveis em tempo real
Para ter uma visão ampla de toda sorte de trabalhos que abordem esse assunto, procurou-se
artigos que reunissem explicativamente diversas técnicas existentes na área, de uma maneira
comparativa e clara. O trabalho Physically based deformable models in computer graphics de
Nealen et al. (2006) faz exatamente isso, relatando e descrevendo os principais trabalhos das
últimas décadas, expondo as peculiariedades e diferenças de cada método de maneira imparcial,
oferecendo ainda, várias fontes bibliográficas bem conceituadas em cada tópico.
2.3.3 Métodos com malha e sem malha
Para efeitos comparativos entre essas abordagens, observou-se trabalhos que houvessem
realizado uma comparação conceitual e qualitativa dos resultados de trabalhos em cada classe.
Isto foi encontrado em To mesh or not to mesh. that is the question de Idelsohn e Onate (2006),
que levanta e procura responder as principais perguntas que envolvem essa grande diferença
entre os métodos de simulação desta área. Também apresentando e descrevendo pontualmente
uma série de trabalhos em cada classe, destacando, por fim, o potencial das técnicas sem malha
que vem sendo desenvolvidas nos últimos anos.
2.3.4 Dos métodos com malha
Para profunda compreensão de todos os fatores que envolvem e caracterizam as técnicas
que utilizam malha, procurou-se os trabalhos mais clássicos e bem estabelecidos na área, em
especial os que tratam sobre método dos elementos finitos, FEM, destacado como principal
modelo. Inicialmente, para entender o que é e como se constrói as malhas tridimensionais
para discretização do contínuo, há o artigo High-quality multi-tissue mesh generation for nite
element analysis de Foteinos e Chrisochoides (2011), que expõe as dificuldades e problemas
em se criar tal estrutura e oferece uma técnica de geração de malhas otimizadas para FEM.
De caráter mais fundamental, a tese On vertex-vertex systems and their use in geometric
and biological modelling de Smith (2006) trata diretamente de triangulação em estruturas de
16
vértices, oferecendo a definição formal de malha.
E para acesso à modelagem e simulação com uso de malha, a obra Introduction to Finite
Element Method de Felippa (2001) que conceitua e expliqua FEM de maneira detalhada e ex-
tensa.
Alternativo ao FEM, a tese Cardiac Surgery Simulation - Graphics Hardware meets Con-
genital Heart Disease de Mosegaard (2006) elabora a modelagem completa de um sistema
massa-mola para simulação cirúrgica, com amplas considerações sobre o modelo computa-
cional e a parte técnica de realizar o processo em GPUs.
2.3.5 Dos métodos livres de malha
Para entender a natureza e motivação dos métodos sem malha, relevou-se os trabalhos
Element-free galerkin methods, Dynamic fracture using element-free galerkin methods e Smoothed
particle hydrodynamics - theory and application to non-spherical stars, respectivamente de Be-
lytschko, Lu e Gu (1994) e Gingold e Monaghan (1977) devido ao pioneirismo e introdução às
características fundamentais dos métodos livres de malha, sobre as quais a maioria dos trabalhos
seguintes se baseou. Demonstrando as possibilidades da então nova abordagem.
Mais recentemente, dá-se grande destaque a Point based animation of elastic, plastic and
melting objects, de Müller et al. (2004), por aplicar efetivamente em um modelo a teoria de
mecânica contínua com os princípios descritos por Belytschko, Lu e Gu (1994) unidos ao SPH
de Gingold e Monaghan (1977), com adendo de uma nova técnica para representação da superfí-
cie, caracterizando mesmo a renderização como livre de malha, através de um sistema baseado
em partículas esparsas que interagem entre si sem qualquer ligação explicitada previamente,
antagônica aos elementos do FEM.
Müller et al. (2005) propõe Meshless deformations based on shape matching, como um
alternativo de caráter mais geométrico do que físico, com uma técnica livre de malha que busca
restaurar a forma inicial definida por pontos em um objeto, focado primariamente na satisfabil-
idade visual, resultando em elevado desempenho, embora descarte realismo.
17
3 Mecânica Contínua
Simular algo requer avaliação e formalização dos fatores que compõe o processo em foco.
Como este trabalho não permeia qualquer assunto inexplorado, parte-se do acesso ao estudo
de Mecânica Contínua, campo da física que engloba teorias, leis e definições em dinâmica de
sólidos deformáveis.
Este capítulo explana conceitos básicos, cuja familiarização é fundamental para posterior
entendimento, crítica e tomada de decisão acerca dos modelos e abstrações computacionais
aplicados neste tipo de simulação.
Contudo, decidiu-se limitar o escopo deste trabalho à dinâmica de materiais que reagem
linearmente às deformações, abstração seguida nas subsequentes definições e modelo computa-
cional aplicado. Nenhum estudo sobre comportamentos não-lineares foi realizado.
3.1 Deslocamento
A ação de forças sobre um objeto com rigidez finita pode modificar seu desenho original.
Portanto, definir sua deformação é o ponto central a ser considerado no modelo.
A princípio, não é conveniente elaborar uma descrição da forma de um sólido e mapeá-la
diretamente a uma nova configuração. Para compor a definição, considerar-se-á o deslocamento
de toda ínfima porção do corpo.
Denominaremos de X as coordenadas de um ponto interior ao objeto B no espaço euclidi-
ano tridimensional, conhecidas como coordenadas de material, ou coordenadas Lagrangianas.
O conjunto destas coordenadas compõe a configuração de referência. (BOWEN, 2010)
Uma deformação causará o reposicionamento destes pontos ao longo do tempo. Trans-
formando as coordenadas de material X para coordenadas espaciais x , também chamadas de
coordenadas Eulerianas. Expresso como função, compondo a equação constitutiva referente
ao deslocamento:
18
x = F(X, t) (3.1)
Sua natureza bijetora relaciona-se com o postulado clássico sobre a matéria não poder ocu-
par mais de um lugar ao mesmo tempo. Portanto, admite-se a existência de sua inversa:
X = F−1(x, t)
A partir disso, calcula-se o deslocamento como um vetor que leva uma partícula da posição
de referência à posição espacial no instante t:
u = F(X, t)−X (3.2)
Essa primeira leitura do estado do objeto permite prosseguir com a análise dos efeitos decor-
rentes da deformação, dando os parâmetros iniciais para moldar as equações constitutivas de um
contínuo. (HEINBOCKEL, 2001)
3.2 Gradiente do Deslocamento
Quando se aplica a análise do deslocamento em um volume não se trata mais de uma sim-
ples diferença de posição, como seria em uma observação unidimensional, pois o deslocamento
passa a ter componentes em todas as direções relativas às três dimensões de espaço. Portanto,
é necessário conhecer as derivadas parciais que definem sua variação em determinado ponto a
ser observado para obter as medidas do deslocamento. Dessa observação surge a definição do
gradiente do campo escalar do deslocamento: a derivação parcial da equação constitutiva 3.1
sobre as três dimensões espaciais.
Fisicamente, o gradiente do deslocamento em um ponto X dá a direção e taxa da maior
variação do deslocamento em relação à vizinhança. Utilizando a seguinte notação:
5u
Análogo ao gradiente de qualquer função, diferencia-se parcialmente o campo escalar do
deslocamento sob cada dimensão de espaço, formando um vetor de três componentes:
5u = {u,v,w} (3.3)
19
Com isto pode-se descobrir a variação do deslocamento ao longo de um vetor arbitrário a
partir deste ponto, apenas projetando o gradiente sobre um vetor~a, como um produto escalar:
5u⊗~a (3.4)
O que resulta em um campo vetorial usado para caracterizar um tensor.
Devido à natureza infinitezimal de um contínuo não há uma função explícita computável
para a obtenção do gradiente do deslocamento (BOWEN, 2010). Resultados de análises numéri-
cas precisam de uma aproximadora que possa ser diferenciada. Isto pode ser feito de maneiras
diferentes de acordo com o modelo aplicado, sendo um dos principais pontos de discriminação
entre os métodos que usam malha e os livres desta. Este trabalho estuda e aplica a elaboração
de uma destas, o método dos mínimos quadrados, ou MLS como o acrônimo do termo em inglês
Moving Least Squares.
3.3 Tensão
A consequência direta do deslocamento de pontos de um corpo é a variação de tamanho
nas diferentes dimensões, o resultado disso caracteriza a tensão como a razão entre as medidas
originais e as deslocadas. Porém, é necessário encaixar este conceito no modelo de descrição
de um contínuo tridimensional.
Ao invés de fazer essa análise pontualmente, como para descrever o deslocamento, pre-
cisamos de uma parcela do corpo como referencial. Comparando o segmento de reta entre dois
pontos A e B num estado não deformado com suas novas posições, assim obtendo a proporção
da deformação nesta direção.
Figure 3.1: Variação da distância entre dois pontos de um contínuo após sofrer deformação.
A figura 3.1 ilustra um corpo sofrendo deformação e a variação da distância entre dois
20
pontos de amostragem do contínuo. Sendo L a distância original definida por AB, e L′ a distância
posterior ao deslocamento. O módulo da tensão nesta orientação é, portanto:
|ε|= L′−LL
(3.5)
Porém, esta análise precisa ser mais elaborada, para que se descubra a tensão com relação às
três dimensões espaciais em determinado ponto. É aqui que usa-se o gradiente do deslocamento.
Tensão passa a ser um tensor com seus componentes definidos por:
εi j =12[ui j +u ji] (3.6)
Com ui j representando os componentes do gradiente do deslocamento. Os quais são obtidos
pelas derivadas parciais do deslocamento sob cada dimensão.
A completa dedução e provas desta definição de tensão podem ser encontradas com detalhes
no livro de Heinbockel (2001).
3.4 Estresse
A lei de Hooke afirma que a deformação de um material elástico causa uma força direta-
mente proporcional ao módulo da tensão. Isso implica que a quantidade de força por unidade de
área em um corpo é relacionada à tensão através de uma função linear, isto é o que caracteriza
o estresse.
Utilizando o gradiente da deformação, numa relação linear com parâmetros obtidos empiri-
camente sobre os materiais, obtém-se o tensor de estresse de acordo com a equação constitutiva.
σ =Cε (3.7)
C é um tensor que descreve a agregação de força a cada componente da tensão.
É conveniente perceber uma característica da tensão, a da simetria de sua matriz.
Observando:
21
ε11 ε12 ε13
ε21 ε22 ε23
ε31 ε32 ε33
(3.8)
Como estas componentes foram obtidas pela diferenciação parcial em cada dimensão, formam-
se as seguintes equivalências:
ε12 ≡ ε21
ε31 ≡ ε13
ε32 ≡ ε23
Justificado pela sua própria definição em 3.6:
12 [ui j +u ji] =
12 [u ji +ui j]
Portanto, as componentes relevantes deste tensor podem ser expressas como:
ε11
ε12
ε13
ε22
ε23
ε33
Para mapear as características físicas às componentes de estresse, incorremos na necessi-
dade de um tensor 6X6, composto por parâmetros que podem ser obtidos empiricamente em
experimentos medindo as propriedades de diferentes materiais.
A relação linear na equação 3.7 é explicitada da seguinte forma:
σ11
σ12
σ13
σ22
σ23
σ33
=
C11 C12 C13 C14 C15 C16
C21 C22 C23 C24 C25 C26
C31 C32 C33 C34 C35 C36
C41 C42 C43 C44 C45 C46
C51 C52 C53 C54 C55 C56
C61 C62 C63 C64 C65 C66
.
ε11
ε12
ε13
ε22
ε23
ε33
(3.9)
Novamente fazendo uso da simetria das componentes direcionais, também existente no
22
tensor de estresse.
No livro de Heinbockel (2001) há extensas deduções acerca de planos de simetria, questão
dependente de cada material. No caso de materiais isométricos, nos quais a tensão se comporta
uniformemente em todas as direções, surge o eixo de simetria:
σ11
σ12
σ13
σ22
σ23
σ33
=
C11 C12 C13 0 0 0
C21 C22 C23 0 0 0
C31 C32 C33 0 0 0
0 0 0 C44 0 0
0 0 0 0 C55 0
0 0 0 0 0 C66
.
ε11
ε12
ε13
ε22
ε23
ε33
(3.10)
Com isso pode-se calcular o efeito de forças sobre um objeto, e obter as resultantes de sua
resistência elástica. Os parâmetros que configuram as propriedades elásticas serão abordados
no capítulo que trata da implementação do simulador, onde será apresentado o tensor sugerido
para a técnica implementada neste trabalho.
23
4 Modelos Computacionais
O capítulo 3 resumiu os principais conceitos para a dinâmica dos sólidos deformáveis.
Agora o modelo matemático descrito deve ser encaixado em uma série de técnicas em ambiente
computacional, balanceando realismo e performance.
O primeiro impasse aparente é exatamente a descrição de um contínuo. Por limitações de
memória e processamento, a análise da deformação não pode ser feita de maneira contínua
para todo o volume de um corpo. Portanto, as equações devem ser aplicadas a uma represen-
tação aproximada dos objetos reais, que reduza para unidades finitas e pontuais os efeitos do
deslocamento ao surgir de uma deformação.
Aqui ocorre a ramificação de técnicas que buscam essa modelagem em duas grandes classes:
com malha e sem malha, doravante referidas pelos termos em inglês mesh e meshless.
4.1 Técnicas Mesh
Mesh é o termo em inglês para malha poligonal, composta por uma coleção de vértices cujas
ligações definem as arestas e sua sequência define os polígonos. Exemplificado pela figura 4.1,
onde os pontos representam os vértices e as linhas suas conexões como arestas. (SMITH, 2006)
Figure 4.1: Malha definida por polígonos com seus vértices e arestas
24
Figure 4.2: Malhas de polígonos representando corpos tridimensionais. (CLKER, )
Em computação gráfica esta é a representação mais típica de objetos, cujas superfícies são
renderizadas como um conjunto de retângulos e triângulos definidos pela malha que envolve
seus volumes, como pode ser visto na figura 4.2 que exibe objetos representados como modelos
de arames.
Estas representações geométricas assemelham-se à discretização de um volume contínuo,
pois uma superfície de complexidade infinita é reduzida a uma aproximação visual passível de
ser computada e manipulada por operações algébricas em seus componentes finitos.
De maneira semelhante, o interior dos corpos pode ser aproximado encaixando poliedros
regulares, assim criando uma malha que preencha o volume com uma estrutura clara de conexões
explícitas. O cálculo do sistema de equações diferenciais parciais passa a ter sua complexidade
reduzida para estas formas, sobre as quais obtém-se a leitura das propriedades físicas de tensão
e estresse decorrentes do deslocamento de seus vértices. Por fim, propagando pelas conexões
as forças resultantes da análise em cada intervalo de tempo.
Com base neste tipo de estrutura se desenvolveram técnicas clássicas para simulação física,
como o Massa-Mola, o método dos elementos finitos, de acrônimo FEM, e suas variações.
Outros exemplos incluem: Método das Diferenças Finitas, Método do Volume Finito, Método
das Fronteiras Finitas. (NEALEN et al., 2006)
Para demonstração de conceito e detalhamento da composição e uso de malha, as duas
primeiras técnicas serão resumidas a seguir.
4.1.1 Massa-Mola
Uma visualização típica do efeito elástico quando se explica lei de Hooke é a de um corpo
preso a uma mola de massa desprezível. Uma força aplicada ao primeiro causa a compressão
25
ou distensão da mola, como ilustrado em 4.3. Assim energia potencial elástica é acumulada,
gerando uma força de resposta diretamente proporcional ao módulo da tensão. Aplicar esta
ideia a uma malha caracteriza o método massa-mola, comparando os vértices com os objetos
dotados de massa e as ligações entre eles atuando como molas ideais. Isso confere um caráter
bastante intuitivo a esta técnica.
L'
L
fi
Figure 4.3: Compressão sofrida por uma mola em decorrência da aplicação de uma força fi,comum deslocamento equivalente a L−L′.
As massas são consideradas como partículas, sobre as quais atuam as forças internas prove-
nientes das conexões por mola com outras massas, como ilustra a figura 4.4. Por isso dá-se
também o nome de Sistema de Partículas a este método, o que pode acabar sendo confundido
com sistemas de partículas sem malha, técnica efetivamente implementada neste trabalho. Por-
tanto, o termo Massa-Mola é o que será utilizado.
26
Figure 4.4: Exemplo ilustrando massas conectadas entre si por molas.
As partículas destes sistemas tem seu movimento regido pela segunda lei de Newton.
fi = mixi
A força atuando sobre uma massa i devido a uma ligação com uma massa j é:
fi = ks(|xi j|− li j)xi j
|xi j|
sendo ks a constante da mola, xi j o vetor da diferença de posições entre i e j, li j o módulo
da distância original, e xi j|xi j| o vetor unitário que dá a orientação da força.
Alguns casos procuram reproduzir comportamento viscoelástico que amortece parte da en-
ergia elástica resultante, isto pode ser feito multiplicando as diferenças de velocidade entre
partículas por um fator de damping, ou amortecimento, mas que deve ser projetado no vetor
que separa as massas a serem comparadas para evitar a perda de velocidade de rotação ou de-
mais efeitos de movimento.
Este é um exemplo de modelo mais canônico para sistemas massa-mola, mas há, ainda,
os que procuram maior generalidade. Alguns destes estabelecem energias de deformação que
tentam ser mantidas para o sistema como um todo, de acordo com algumas restrições. Nestes
as forças são consideradas como derivadas destas energias sobre as posições de cada partícula.
Em resumo, o massa-mola requer a resolução de um sistema de equações diferenciais or-
27
dinárias, o que já é contemplado pelo esquema de integração numérica, etapa presente em qual-
quer simulação física, e que será detalhada posteriormente.
Aplicações práticas deste método comumente abordam simulação de tecido, tentando repro-
duzir o efeito elástico causado pela estrutura de fibras, algo ao qual a ideia de malha composta
por massas e molas encaixa-se bem, geralmente com resultados visuais bastante satisfatórios.
Figure 4.5: Simulação em tempo real de tecido usando massa-mola.
A simplicidade de implementação dessa técnica unida à sua eficiência computacional lhe
confere um bom apelo para aplicações visuais que não dependem muito de acurácia, pois este
método possui a desvantagem de não se basear na teoria da elasticidade em sua totalidade, con-
templada pelas definições no capítulo 3. Apesar de usar a lei de Hooke, o comportamento elás-
tico de um contínuo não se resume a este tipo de transmissão de forças. Isto o torna altamente
sensível à granularidade e topologia da malha e seu resultado não converge para o comporta-
mento real quando se aumenta a complexidade da estrutura. Já evidenciando empecilhos desta
classe de métodos dependentes de malha.(NEALEN et al., 2006)
Apesar das limitações desse método, a tese de Mosegaard (2006) apresenta a implemen-
tação em GPU de simulação de cirurgia cardiaca com massa-mola com ótimos resultados vi-
suais e de elevada eficiência.
28
4.1.2 FEM - Método dos Elementos Finitos
FEM consolida-se como um dos métodos mais populares para se resolver equações difer-
enciais parciais, (NEALEN et al., 2006), aqui referenciadas como o acrônimo PDE do inglês
Partial Diferential Equation, usado para discretizar volumes contínuos e permitir a computabil-
idade para várias aplicações, em especial a simulação de corpos deformáveis.
Ele prossegue com a definição de uma malha tridimensional irregular, cujos poliedros car-
acterizam os mencionados elementos finitos, ilustrado em duas dimensões na figura 4.6. E sobre
eles se aplica a equação diferencial de elasticidade:
ρ x = ∇σ + f
,
sendo ρ a densidade do material e f a resultante das forças externas.
Os cálculos passam a ser realizados sobre os nodos dos elementos, avaliando a variação de
suas posições x(t). Mas para obter o gradiente precisamos aproximar a função x(m, t) descritora
do contínuo, usando os valores dos nós:
x(m, t)'∑xi(t)bi(m)
,
onde bi() são funções base fixas que valem 1 no nodo i e 0 nos demais, conhecida como a
propriedade Delta de Kronecker, e m a massa. (NEALEN et al., 2006)
Figure 4.6: Elemento 2D sofrendo deformação e as forças elásticas atuando nos nodos.
Através da função aproximada, o espaço infinito de soluções é reduzido a um subespaço
finito, assim correlacionando a discretização do volume à essa sub imagem das PDE originais.
Contudo, nenhuma função aproximadora é capaz de resolver as PDE originais com exatidão,
29
pois haverá erro residual decorrente do afastamento da solução aproximada à ideal.
No método de Galerkin é realizada a busca pela aproximação que minimize esses resí-
duos, de forma que extraia-se a melhor precisão alcançável pelo modelo discreto do volume
(HUNTER; PULLAN, 2001). Este tipo de aproximação será detalhado nos capítulos que abor-
dam a implementação deste trabalho.
À parte das considerações básicas sobre FEM, existem implementações distintas. Aqui se
resume o FEM Explícito, uma das abordagens mais simplistas.
No FEM Explícito as massas e forças externas são concentradas nos vértices, o que acaba
se assemelhando ao massa-mola, pois os nodos atuam como massas e suas conexões com os
nodos adjacentes passam a representar molas. Quando um elemento sofre deformação as forças
atuantes nos seus nodos podem ser calculadas com a função aproximadora, obtendo-se tensão e
estresse a partir do campo de deslocamento. Então a energia de deformação de um elemento é
dada por:
E =
ˆε(m) ·σ(m)dm
As forças podem então ser calculadas como a derivada da energia em relação às posições
dos nodos. E linearizando a relação das forças conectadas a um elemento e fica como:
fe = Keue
,
onde Ke é a matriz de rigidez do elemento e como as forças vindas de elementos adja-
centes se aglomeram pode-se formar uma matriz de rigidez somando das matrizes de todos os
elementos:
K = ∑Ke
.
K terá dimensão de 3n× 3n, e para expandir as matrizes individuais a esta dimensão
preenche-se com zeros nas posições referentes a nodos não adjacentes ao elemento. Finalmente,
usando a elasticidade linearizada a equação do movimento para toda a malha se torna:
Mu+Du+Ku = fext
30
,
onde M é a matriz de massas, D é a matriz de amortecimento, ou damping, e fext as forças
externas.
Uma introdução completa do método dos elementos finitos encontra-se no livro de Felippa
(2001). E para demonstração de aplicabilidade direta o trabalho de Zachow et al. (2000) apre-
senta a utilidade do FEM para simulação cirúrgica em seu caso de estudo de planejamento de
reconstrução facial assistido por simulação física da mesma natureza aqui estudada.
Diferente do massa-mola, o FEM, por seguir a teoria da elasticidade, tende gradualmente
ao comportamento real conforme se refina a discretização rumo ao contínuo ideal, sacrificando
eficiência ao longo disto. Aplicações científicas e de engenharia, como para análise estrutural,
costumam fazer uso disso, abrindo mão de resultados em tempo real para agregar precisão, mas
de maneira ainda viável, executando FEM sobre malhas de granularidade fina.
A realização prática deste refinamento introduz os principais desafios e desvantagens do
uso de FEM, devido a sua sensibilidade à estrutura da malha e a grande dificuldade de produzir
e modificar a mesma. O que motiva a busca por uma alternativa que iguale ou supere suas
qualidades e livre-se dos problemas inerentes à geração de manutenção da malha.
4.1.3 Construção e preservação da Malha
O principal desafio enfrentado por esta classe de métodos, em especial o FEM, parte da ne-
cessidade primordial de construir-se uma estrutura de malha consistente, pelas peculiariedades
dos casos a serem simulados, e do método utilizado.
Situações diferentes costumam requerer modelagens específicas, especialmente em casos
delicados como o uso de massa-mola, cujas ligações entre as massas altera drasticamente o
resultado. Isto indica que sua construção está longe de ser trivial, sendo a malha por definição
uma série de ligações explícitas entre nodos. Por exemplo, a preocupação da distruibuição das
forças em grades de massas e molas, seja para dar uniformidade ou reforçar alguma direção,
como mostra a estrutura ilustrada na figura 4.7 que define todas as conexões possíveis entre os
vértices para resultar em comportamento isotrópico, que é a tensão com mesmo comportamento
em todas as direções.
31
Figure 4.7: Malha em sistema de massa-mola planejada para comportamento isotrópico
(NEALEN et al., 2006)
Mesmo usando FEM não basta preocupar-se apenas com a subdivisão do espaço, pois as
características dos poliedros resultantes influenciam a estabilidade numérica ao longo de todo o
processo. Portanto, além da preocupação em encaixar poliedros diversos em formas tridimen-
sionais arbitrárias, isto também não pode ser feito de maneira ingênua, pois além da ligação
entre os elementos, formar poliedros com ângulos muito pequenos ou muito grandes causa er-
ros de interpolação e mau condicionamento das matrizes de rigidez.
Existem trabalhos específicos sobre a geração de malhas otimizadas para FEM, como o
de Foteinos e Chrisochoides (2011), que preocupa-se em evitar a formação de poliedros com
ângulos inadequados, eliminando-os em troca da adição de seu circuncentro, centro da menor
esfera que o engloba, e com a fidelidade da superfície com técnica de rejeição de tais pontos
adicionados de forma a excetuar os gerados próximos à superfície. Malhas resultantes desta
técnica podem ser vistas nas figuras 4.8 e 4.9.
Não bastando as considerações e esforços aplicados na geração da malha, o método de sim-
ulação utilizado deve preocupar-se também em manter a coerência das ligações entre os nodos
após sofrer sucessivas ou grandes deformações, o que, dependendo do caso, pode adicionar
custo de processamento à cada iteração da simulação.
Este tipo de obstáculo imposto pela necessidade de uma malha estimulou a pesquisa de
métodos que não necessitem de tal estrutura, motivando o interesse deste trabalho.
4.2 Técnicas Meshless
Pesquisas que buscam respostas para contornar os empecilhos inerentes ao uso de malha
começaram a surgir, especialmente pelo desejo de obter um interpolante polinomial através de
32
Figure 4.8: Malha representando estruturas anatômicas de um tronco humano, com a dis-tribuição dos ângulos dos poliedros. (FOTEINOS; CHRISOCHOIDES, 2011)
Figure 4.9: a) Malha com buracos. b) Fidelidade da superfície aprimorada com técnica derejeição de pontos. c) Corte transversal mostrando a estrutura interna da malha. (FOTEINOS;CHRISOCHOIDES, 2011)
33
pontos de amostra que minimize a distância aos pontos desconhecidos do contínuo. O trabalho
pioneiro nesta busca foi o de Nayroles, Touzot e Villon (1992), que trouxe o desenvolvimento
de uma aproximação mais suave e contínua que necessitasse apenas de pontos de controle e não
de ligações explícitas entre elementos como poliedros. Suas ideias procuravam estender o FEM
para um caráter mais geral e lançaram a base para esta linha de pesquisa.
Poucos anos depois, uma aplicação direta desse tipo de técnica foi proposta por Belytschko,
Lu e Gu (1994) com o chamado Element-free Galerkin Method, EFGM. Métodos de Galerkin,
na matemática, referem-se à conversão de problemas com operadores contínuos em discretos,
como a já mencionada aproximação do contínuo pelo método FEM. O EFGM faz o mesmo,
porém, com o uso de um interpolador chamado Moving-Least-Squares, método dos mínimos
quadrados, ou apenas MLS. Este trabalho apresenta um capítulo que descreve especificamente
este interpolador. Com o exemplo da figura 4.10 do trabalho de Belytschko, Lu e Gu (1994),
onde a simulação de rachadura em um plano precisou apenas de uma grade de pontos sem
qualquer conexão explícita.
Figure 4.10: Simulação de rachadura com EFGM, onde o domínio é representado apenas por
pontos, sem necessidade de ligações. (BELYTSCHKO; TABBARA, 1996)
Não apenas acerca da interpolação, houve também uma mudança na consideração da com-
posiçao das unidades elementares influenciado diretamente pelo trabalho de Gingold e Mon-
aghan (1977) com o método Smoothed particle hydrodynamics, ou SPH, muito usado para
dinâmica de fluidos onde faz-se o uso da representação das partículas como componentes de
um contínuo, exemplificado pela figura 4.11 que retrata momentos na simulação de um fluido
viscoso com o protótipo resultante do trabalho de NOBREGA (2010). Extensões desta ideia
também foram aplicadas à análises estruturais, formando a base para vários dos métodos livres
de malha. Pois a representação de partículas encaixa-se nas características dos interpoladores
34
suaves usado por esta classe de métodos. (NEALEN et al., 2006)
Figure 4.11: Dinâmica das partículas em protótipo de simulação de fluídos baseado em SPH.
(NOBREGA, 2010)
A partir destes trabalhos surgiram diversas linhas de pesquisa que acreditam no potencial
dos métodos sem malha, inclusive com a demonstração de sua funcionalidade em animações
com alto grau de realismo, como o trabalho de Müller et al. (2004) que apresenta uma técnica
para animação elástica e plástica inteiramente baseada nos princípios de mecânica contínua
aplicado a abstrações de partículas, herdando os princípios do SPH unidos à base matemática
do EFGM de interpolação polinomial para aproximação do contínuo. Compondo, então, um
método cuja diminuição da granularidade tenda ao real.
Destaca-se, pelos resultados diretos, sua técnica de Shape-matching, Müller et al. (2005),
focada em agradabilidade visual com elevada eficiência computacional, mas que diverge em
realismo das técnicas que apliquem as teorias físicas, pois utiliza de uma ideia geométrica de
restauração da forma inicial, sendo capaz de suportar grandes deformações mantendo estabili-
dade incondicionalmente.
Decidiu-se a implementação da técnica baseada em pontos para avaliar sua eficiência e es-
tabilidade rumo a aplicações em tempo real passíveis de serem utilizadas em áreas científicas,
em especial a simulação de tecidos e orgãos biológicos. Portanto, os capítulos seguintes abor-
dam todos os princípios e detalhes pertinentes a esta técnica, em paralelo com a implementação
da mesma em um protótipo de simulador interativo de corpos deformáveis.
35
5 Discretização do contínuo
Quando se trata de um problema de movimento no espaço, como na cinemática, ou de in-
teração entre forças, como na mecânica clássica, os elementos envolvidos podem ser compreen-
didos como partículas, cujos comportamentos são descritos analiticamente. Assim, o domínio1
de tais definições abrange todo o universo estudado. Ou seja, aplica-se análise absoluta direta-
mente sobre a partícula, indiscriminando sua posição no espaço e tempo. Na matemática, este
tipo de análise classifica-se como uma formulação forte.
Um modelo analítico dessa natureza é inviável na dinâmica dos sólidos, pois o contínuo
deve ser considerado como um todo, ao longo do qual a análise diferencia-se. O que indica a
necessidade de uma solução numérica que parta do princípio de integrar subdivisões do domínio
(BOWEN, 2010).
Nas equações da dinâmica dos sólidos, necessita-se de uma função diferenciável do deslo-
camento para poder obter o gradiente deste, então as técnicas livres de malha buscaram encon-
trar uma aproximação válida para essa função ideal.
5.1 Método de Galerkin
No método dos elementos finitos o pioneirismo de Boris G. Galerkin é aplicado para a
subdivisão do espaço, onde a função do campo contínuo do volume é trocada por funções lin-
eares de base, que definem a aproximação paramétrica dos elementos finitos. Cada subdivisão
é chamada de elemento, e a soma das integrais de cada elemento aproxima-se da integral do
contínuo, porém com um erro residual.
A suposição de Galerkin força com que o erro residual seja ortogonal ao espaço das funções
base, o que faz com que o erro seja reduzido monotonicamente conforme a malha é refinada.
Ou seja, quanto maior o número de elementos dividindo o espaço maior a semelhança com o
preenchimento contínuo, e colateralmente a complexidade da simulação (HUNTER; PULLAN,
1aqui o termo “domínio” faz analogia direta com o domínio de funções
36
2001).
O início das pesquisas em modelos livres de malha aproveitou o conceito dos métodos
de Galerkin (BELYTSCHKO; TABBARA, 1996), almejando resultados com as mesmas car-
acterísticas, em especial sobre o comportamento do erro residual. Porém, livrando-se do em-
pecilho na criação dos elementos geométricos e sua conectividade que compõe a malha. As
pesquisas começaram a usar formulações que se aproximassem do contínuo com base em um
conjunto de amostragem, geralmente pontos. E a técnica encontrada em comum na maioria
desses métodos é a dos mínimos quadrados, especificamente Moving Least Squares (NEALEN
et al., 2006).
5.2 Amostragem e função aproximadora
A aplicação da ideia dos métodos de Galerkin na dinâmica dos sólidos consiste em substituir
o espaço contínuo de um volume por uma aproximação composta pela integralização de porções
do corpo. O que no FEM, principal representante das técnicas com malha, é feito com poliedros
interconectados (NAYROLES; TOUZOT; VILLON, 1992).
Nos métodos livres de malha, a preocupação em contornar os problemas de criação e
manutenção da mesma não exime da necessidade de obter unidades de amostragem de um
corpo. Esta passa a consistir de um conjunto N de pontos arbitrariamente distribuídos no
domínio do volume, sem qualquer relação de conexão explícita que não a de vizinhança. Por-
tanto, a geração desta amostra consiste em registrar as posições xi no espaço R3de cada ponto
do conjunto.
O próximo passo é compor uma função f (x) de abrangência global que se aproxime dos
valores escalares fi nos pontos de N. Sendo desejável as propriedades de ortogonalidade entre
o erro residual e o espaço da função. Para alcançar este requisito, f (x) deve ser obtida por meio
da minimização do erro entre a aproximação e o valor da amostra f (xi)− fi . Existe uma série
de técnicas familiares que buscam a minimização da soma dos quadrados deste erro, explicadas
a seguir.
5.3 Least Squares
Esta é a técnica que serve de base para a minimização do somatório do erro E definido por:
37
E = ∑i‖ f (xi)− fi ‖2 (5.1)
Onde queremos encontrar a função f pertence ao espaço ∏dm de polinômios de grau m e
dimensão d, que tenha menor erro residual possível. Ela pode ser descrita como:
f (x) = b(x)T c (5.2)
Sendo b(x) = [b1(x), ...,bk(x)]T o vetor da base polinomial que define o tipo de aproxi-
madora que deseja-se construir em cada aplicação, formando um polinômio de qualquer grau
e dimensões necessários. Como exemplo, para gerar uma função que descreva em forma de
superfície quadrática a aproximação de uma série de pontos, usa-se a base de segundo grau e
duas dimensões b(x) = [1,x,y,x²,xy,y²]. Ou, como no contexto deste trabalho, para modelar
a função aproximadora linear do deslocamento de pontos usando base de primeiro grau e três
dimensões b(x) = [1,x,y,z].
Com a base definida, resta encontrar o vetor de coeficientes c = [c1, ...,ck]T do polinômio
aproximador que minimize a equação 5.1. Isto consiste diretamente em encontrar o mínimo de
uma função, o que é feito procurando os valores que levem as derivadas parciais a zero.
Usando a equação 5.2 em 5.1, temos cada valor em c como uma incógnita, sobre cada qual
diferencia-se parcialmente compondo ∇E = [∂E/∂c1, ...,∂E/∂ck] . Igualando cada termo do
vetor ∇E a zero, forma-se um sistema linear de coordenadas. As derivadas parciais de E sobre
cada componente c se apresentam da seguinte forma:
∂E/∂c1 = ∑i
2b1[b(xi)T c− fi] = 0
∂E/∂c2 = ∑i
2b2[b(xi)T c− fi] = 0
...
∂E/∂ck = ∑i
2bk[b(xi)T c− fi] = 0
Que pode ser escrito em notação de matriz:
∑i
2b(xi)[b(xi)T c− fi] = 2∑
ib(xi)[b(xi)
T c− fi] = 0
38
Agora dividindo pela constante em ambos os lados e distribuindo a multiplicação:
∑i[b(xi)b(xi)
T c−b(xi) fi] = 0
E, finalmente, isolando c:
∑i
b(xi)b(xi)T c = ∑
ib(xi) fi
c = [∑i
b(xi)b(xi)T ]−1
∑i
b(xi) fi (5.3)
Se o determinante da matriz ∑i b(xi)b(xi)T for zero, significa que ela é singular e não
inversível. Portanto incapaz de ser utilizada para formar a equação 5.2. Do contrário, aplicando
os coeficientes 5.3 no polinômio 5.2 forma-se uma função aproximadora para o domínio dos
pontos amostrados através dos mínimos quadrados (NEALEN; DARMSTADT, 2004).
Figura 5.1: Dois exemplos de superfícies geradas a partir de duas amostras de pontos aproxi-madas pelos polinômios gerados por mínimos quadrados. (NEALEN; DARMSTADT, 2004)
A justificação visual da aplicação dessa técnica pode ser vista na figura 5.1, onde a partir
de amostras relativamente escassas pode-se gerar superfícies que representam a informação
contínua que se aproxima dos pontos do domínio, usando mínimos quadrados para encaixá-los
em polinômios quadráticos bidimensionais.
39
5.4 Weighted Least Squares
A aproximadora anteriormente definida tem abrangência global no domínio da amostra, e
considera a influência de todos os pontos igualmente. Para casos em que se deseja diferenciar
uma função do domínio de maneira a considerar a relevância de cada ponto, a técnica dos
mínimos quadrados foi adaptada para incluir uma função de peso (NEALEN; DARMSTADT,
2004). O método resultante é o dos mínimos quadrados ponderados, aqui referido pela sigla do
termo em inglês WLS.
Nele uma função de peso w é inserida na minimização do erro:
E = ∑i
w(‖ x−xi ‖) ‖ f (xi)− fi ‖2
Ela atua sobre o valor do erro em relação à distância entre um ponto x e os pontos da
amostra, além de considerar uma região de abrangência r em torno de x. Assim a função
aproximadora gerada deixa de ser válida no sentido global, e passa a ser em relação ao ponto x
que atua como centro local. Os coeficientes do polinômio passam a ser dependentes deste ponto
de análise:
f (x) = b(x)T c(x)
O restante do processo de minimização é análogo ao LS convencional, porém, deve-se gerar
uma aproximadora para cada ponto de análise desejado (NEALEN; DARMSTADT, 2004).
O critério de definição da função w depende de cada aplicação, apenas que ela seja contínua
no domínio r em torno de x. Sua influência atua diretamente no resultado de cada aproximadora
gerada, pois agora os pontos não são considerados de maneira equivalente, o que faz com que
o polinômio seja curvado de acordo com a distância dos pontos vizinhos, como exemplificado
na figura 5.2, onde a influência do ponto A sobre a curva é menor do que a do ponto B que se
encontra mais próximo.
5.4.1 Função de peso
As funções de peso são responsáveis por distribuir a influência em torno do ponto de suporte
de acordo com cada aplicação desejada. Ela pode ser modelada para cada caso, e basicamente
dá a ideia do decaimento da interferência dos pontos conforme sua distância do suporte. Nor-
malmente definida como uma função de duas variáveis:
40
X
A
B
Figura 5.2: A curva em azul representa uma aproximação por mínimos quadrados, onde ainfluência dos pontos é equivalente e global, o polinômio gerado tende a se aproximar de todoscom o menor erro possível. A curva rosada representa uma aproximadora gerada por mínimosquadrados ponderados, onde a distância dos pontos ao centro de análise é considerada. Destaforma a influência de pontos mais afastados sobre a curva é inferior à dos pontos próximos a x.
w(d,r)
Onde d é a distância de cada amostra ao ponto de suporte, e r o raio dessa partição. Para
garantir diferenciabilidade da função aproximadora, w deve ser contínua enquanto d ≤ r, no
geral com w(d,r) = 0 para d > r.
Dois interesses são confrontados na composição desta função: da necessidade de um inter-
polante ou de um aproximador. Se limd→0 w(d,r) = ∞, a função gerada toma características
de um interpolador. E quando w(d,r) decair rapidamente quando d→ ∞, obtem-se um aproxi-
mador local (LEVIN, 1998).
41
Figura 5.3: Em gráfico o comportamento de três diferentes funções de peso com raios de
tamanho 1. A gaussiana em roxo (NEALEN; DARMSTADT, 2004), poly6 em verde (MÜLLER
et al., 2004) e spiky kernel em amarelo (NOBREGA, 2010).
Neste trabalho adotou-se a função poly6, indicada por Müller et al. (2004):
wpoly6(d,r) =
315
64πr9
(r2−d2)3 d < r
0 d ≥ r(5.4)
5.4.2 Partição do domínio
O resultado do WLS é uma função válida localmente, mas pode-se compor uma função de
abrangência global. Isto é feito garantindo que para qualquer ponto no domínio exista uma
função aproximadora que o contenha. Ou seja, deve-se gerar funções dispersas por todo o
contínuo, de forma que a união de suas localidades contemple todo o domínio.
Inserindo esta questão no problema de modelar o volume de um corpo, deve-se espalhar
por todo o domínio pontos de suporte para o WLS. Então a união do domínio das funções de
peso deve resultar, idealmente, no domínio global do corpo. A figura 5.4 demonstra em duas
dimensões uma possível distribuição de pontos de suporte com os domínios de suas funções de
peso cobrindo a maior parte de um objeto. Para a simulação de um volume sem usar malha,
extrapola-se essa ideia para três dimensões, tentando preencher a totalidade do corpo.
Na partição de domínio, usa-se a união dos domínios locais das funções de peso:
Wj =w j(x)
∑k wk(x)
42
Representando a proporção do peso em relação à soma de todos os pesos. Com isso, a
composição de uma aproximadora local geraria:
f (x) = ∑j
Wj (x)b(x)T c(xj)
Em decorrência da soma de todas as partições do domínio por WLS, unindo todos os aprox-
imadores pela integração de suas regiões de influência (NEALEN; DARMSTADT, 2004).
Figura 5.4: Demonstração de partição de domínio com pontos de análise com abrangência de
raios grandes, buscando cobrir toda a área do domínio.
5.5 Moving Least Squares
O método de aproximação Moving Least Squares, aqui referenciado pela sigla MLS, foi
proposto por Lancaster e Salkauskas (1981), para gerar superfícies implícitas a partir de nu-
vens de pontos. Eles realçaram a natureza dupla de interpolação e aproximação conforme a
configuração da função de peso.
Sua ideia representa um fino contraste com a partição do domínio. Ao invés da integral-
ização dos domínios da partição, gera-se uma função com o WLS em um ponto arbitrário do
volume, e move-se para outra posição, gerando funções por todo o domínio, por isso sua denom-
inação de mínimos quadrados móveis. Computacionalmente falando, isso permite o controle do
43
refinamento da simulação. Inclusive de sua natureza, pois alterando a continuidade da função
de peso alterna-se entre interpolação e aproximação (NEALEN; DARMSTADT, 2004).
A técnica de simulação de corpos deformáveis baseada em pontos de Müller et al. (2004)
faz uso deste método, assim como o trabalho pioneiro de BELYTSCHKO e TABBARA (1996),
inspirados pelo bem sucedido trabalho de Lancaster e Salkauskas (1981). Dado o fato de que
utilizam como amostra nuvens de pontos, apesar da diferente aplicação.
44
6 Simulação elástica e plástica baseadaem pontos
Dos frutos da pesquisa em simulação de corpos deformáveis sem malha das últimas dé-
cadas, optou-se pela técnica baseada em pontos de Müller et al. (2004), por seu fundamento
ser íntimo da modelagem física de mecânica contínua. Os autores fornecem demonstrações
de seu funcionamento bem sucedido, incluindo um protótipo executável, dentre uma série de
vídeos e imagens. A intenção é seguir em uma linha de pesquisa que aproxime realismo de
desempenho, pois o interesse está não apenas em criar animações visualmente agradáveis, mas
convergir gradualmente ao realismo de métodos consagrados da classe antagônica baseada em
malha, como o FEM.
Este capítulo explica a idéia do método e detalha todos os passos necessários para realizar
este tipo de simulação. O texto e fórmulas seguem a notação utilizada pelos autores, explicando
devidamente cada termo.
6.1 Visão geral do método
A principal característica deste método é a unificação de duas ideias de simulação sem
malha: as porções do corpo são interpretadas como partículas semelhante às do SPH e cada
partícula é um ponto de suporte na análise com MLS como no EFGM, ambos descritos no
capítulo 4.
O interior de um corpo é representado por um conjunto arbitrariamente distribuído de
partículas denominadas phyxels. Cada phyxel possui um valor de massa, densidade, veloci-
dade, aceleração, raio e uma lista phyxels vizinhos que se encontram na abrangência do raio.
O exterior, ou a superfície a ser renderizada, é representado por pontos chamados de surfels,
que são dotados de posição, vetor normal à superfície, raio e uma lista de phyxels vizinhos. O
método se baseia fortemente em técnicas de computação gráfica com pontos para renderização
45
de superfícies complexas representadas por surfels. Mas para a simulação, a essência é a inter-
ação entre os pontos da superfície e as partículas do interior. Abstrai-se, portanto, dos métodos
de renderização de tais superfícies (GROSS; PFISTER, 2007).
Apesar de ser um método que não utiliza malha com ligações explícitas, existe uma relação
de vizinhança entre os phyxels. A representação do corpo inicia com uma distribuição qualquer
de phyxels pelo volume que deseja-se simular, sendo cada um configurado de acordo com N
vizinhos mais próximos, registrando seus vizinhos para as etapas posteriores da simulação.
Uma vez configuradas as partículas e suas vizinhanças, passa-se a iterar no tempo sobre a
simulação realizando os seguintes passos sobre cada phyxel do corpo:
ut︸︷︷︸deslocamentos
→ ∇ut︸︷︷︸gradiente
→ εt︸︷︷︸tensao
→ σt︸︷︷︸estresse
→ ft︸︷︷︸f orcas
→ ut +4t︸ ︷︷ ︸integracao
Onde cada etapa aplica os cálculos com base na aproximação por WLS gerada tendo o
phyxel corrente como ponto de suporte. A junção do resultado da iteração sobre cada phyxel
resulta na aplicação de MLS, onde a marcha pelo domínio passa por cada partícula gerando um
aproximador por mínimos quadrados ponderados tendo seus vizinhos como amostra.
Ao final de cada iteração, a integração temporal atualiza a posição de cada phyxel repro-
duzindo seu comportamento de partícula, e os surfels tem sua posição atualizada com base nos
phyxels mais próximos. O resultado da sobreposição temporal das iterações é a aproximação
da dinâmica de um objeto deformável com comportamento modelado pelas leis da mecânica
contínua.
As seções a seguir detalham como o método aplica a análise vista no capítulo 3 sobre esta
representação de corpo contínuo em forma de phyxels e surfels, em uma aplicação que simule
objetos deformáveis com resultados em tempo real, passíveis de iteratividade e com grande
apelo visual como pode ser visto na figura 6.1.
Figura 6.1: Demonstração do resultado da simulação de plasticidade e elasticidade com a téc-
nica baseada em pontos, que podem ser vistos na figura da esquerda. (MÜLLER et al., 2004)
46
6.2 Inicialização dos phyxels
Antes de se poder iterar sobre as partículas elas devem ser criadas e configuradas. Tudo
parte de uma distribuição inicial irregular de um número qualquer de phyxels no volume do
corpo. Aqui a quantidade de partículas utilizadas interfere diretamente na complexidade e refi-
namento da simulação, pelas razões já comentadas no capítulo 5.
Em dada distribuição, um phyxel irá corresponder a uma porção do corpo, isso significa
que ele representa uma massa e um volume. Para definir isto, deve-se avaliar sua vizinhança,
pois phyxels com vizinhos mais distantes deverão representar uma área maior e, de acordo com
a densidade desejada, uma massa maior, já que existem poucos outros pontos de amostra em
torno do local. Para gerar essa forma de equilíbrio da representação, avalia-se as partículas
como descrito a seguir.
• Para cada phyxel pi faça:
1. Obtenha distância média d aos k vizinhos mais próximos.
2. Defina o raio ri de pi como um múltiplo de d .
3. Defina a massa m de pi como mi = sd3ρ .
4. Adicione na lista de vizinhança de pi todos os phyxels com distância menor do que
ri.
Os autores sugerem k = 10 e ri = 3d, e na definição da massa ρ é a densidade do material e s é
um fator de escala que aproxime a densidade de todos os phyxels.
6.3 Calculando o gradiente do deslocamento
Como já definido, o gradiente do deslocamento é a derivada parcial do mesmo sobre cada
componente. Mas para realizar isso, precisa-se de uma função diferenciável, quem fornece tal
função é o método MLS , que aplica WLS sobre cada phyxel.
Primeiro é necessário entender como enquadrar o MLS no modelo.
6.3.1 Aplicando MLS
Deseja-se obter as derivadas da função do deslocamento para cada phyxel, mas essa função
não é explicitamente conhecida. Tentando descobrir essa função para cada phyxel usando WLS,
47
a minimização do erro forma:
E = ∑j
(u j−u j
)2 wi j
Com w sendo a função de peso poly6 vista em 5.4.
A minimização deste erro por WLS conduz às derivadas parciais de ui j sobre os coeficientes
da aproximação. Como ui j é uma função sobre as coordenadas de material, que são: a distância
inicial entre o suporte i e cada vizinho j, a base polinomial b consiste em xi j, e os coeficientes
c a serem minimizados correspondem às componentes que levam xi j à nova posição deslocada.
Portanto, o gradiente do deslocamento é a minimização desses coeficientes em si. Fazendo o
comparativo entre as equações:
∑i
b(xi)b(xi)T c = ∑
ib(xi) fi
↓(∑
jxi jxT
i jwi j
)∇u = ∑
jui jxi jwi j
Para obter o gradiente resta isolar:
∇u =
(∑
jxi jxT
i jwi j
)−1
∑j
ui jxi jwi j (6.1)
Para os casos em que o número de phyxels na vizinhança for menor do que 4, incluindo o
próprio, ou eles forem colineares ou coplanares, a matriz A = ∑ j xi jxTi jwi j se torna singular e
não pode ser invertida. Quando este for o caso, ou mesmo quando ela for má condicionada, o
que levaria a componentes grandes demais, deve-se aplicar a inversão através da decomposição
por SVD, descrita em 6.3.2. Ademais, seu cálculo só precisa ser efetuado uma vez ao longo da
simulação, pois ela depende apenas das coordenadas de materiais, que são as posições iniciais
dos phyxeis no corpo. Isto representa uma grande otimização no custo do processo de cada
iteração.
48
6.3.2 Inversão por SVD
Em matrizes singulares, ou cujo mal condicionamento as aproxima da singularidade, tipi-
camente não se encontra uma solução para sua inversa. Porém, negar a inversibilidade de uma
matriz pode significar o descarte de informação útil que está contida em dimensões bem com-
portadas da mesma. Uma maneira de verificar o estado da matriz, que seria perceber os motivos
de sua singularidade, é usando a técnica de Single Value Decomposition (PRESS et al., 1992).
Um teorema da álgebra linear afirma que toda matriz AM×N com número de linhas maior
que o número de colunas pode ser escrita como o produto de uma matriz UM×N de colunas
ortogonais, uma matriz diagonal WN×N e uma matriz transposta VTN×N :
A = U ·W ·VT
U e V são ortogonais, pois suas colunas são ortonormais, isso significa que :
U ·UT = VT ·V = 1
Inverter usando a decomposição por SVD no caso de matrizes quadradas, como a matriz A
gerada pelo MLS, se torna trivial e acontece da seguinte forma:
A−1 = V ·W−1 ·UT (6.2)
Pelas suas propriedades ortonormais, as matrizes U e V tem a inversa dada pela sua trans-
posta, e a matriz W por ser diagonal tem sua inversa como a diagonal com o inverso multi-
plicativo de seus elementos. E é neste ponto que acontece a análise da singularidade ou mau
condicionamento, mas antes precisa-se entender como formar cada uma das três matrizes.
Por A ser quadrada, sua decomposição consiste em separar os auto-valores dos auto-vetores,
que podem ser interpretados respectivamente como a magnitude em cada dimensão, e os vetores
que definem as direções. U e V são compostos por auto-vetores de A ·AT e a matriz W pelos
auto-valores dispostos na diagonal.
Para avaliar o condicionamento ou a singularidade de A, precisa-se observar os valores de
W:
49
W =
w1 0 0
0 w2 0
0 0 w3
Como sua inversa é dada pelo inverso multiplicativo, se um dos valores da diagonal de W
for igual a zero ou muito pequeno, um ou mais termos de W−1 representará uma indefinição,
ou um número muito grande, pois:
limwk→0
1wk
= ∞
O que impediria sua inversão. Mas isto pode ser tratado de forma a preservar o restante da
informação válida. Por exemplo, no caso de phyxels coplanares, pode-se descartar a dimensão
onde ocorre a singularidade. Isto é feito zerando o inverso multiplicativo do termo causador
da indefinição, fazendo 1wk
= 0. O mesmo pode ser feito no caso de phyxels muito próximos
de serem coplanares, em que o inverso resultaria em um número muito grande que causaria
instabilidade numérica na integração temporal. Para isso estipula-se um limite inferior ε para
os auto-valores.
• Para cada auto-valor wk em W:
– Se wk < ε defina 1wk
= 0
Desta forma descarta-se os pontos de singularidade e preserva-se para uso na simulação as
dimensões válidas da matriz A. Como demonstrado na figura 6.2 onde um conjunto de pontos
encontra-se numa distribuição coplanar e exibe-se os eixos das dimensões de tensão.
50
Figura 6.2: Distribuição coplanar de partículas.
Os eixos azul e verde pertencentes ao plano consistem em dimensões válidas com auto-
valores distantes de zero. Já o eixo em vermelho é ortogonal ao plano, assim representando a
dimensão a ser desconsiderada pelo inversão com SVD, por gerar auto-valor igual ou próximo
de zero.
6.4 Calculando a tensão
O método de MLS forneceu o gradiente do deslocamento no phyxel que está sendo anal-
isado. A partir dele a técnica realiza a simulação de duas propriedades dos corpos deformáveis:
elasticidade e plasticidade.
O gradiente do deslocamento conduz ao tensor elástico, e sobre este, de acordo com pro-
priedades do material, o tensor plástico. O cálculo indivual de cada um é visto a seguir.
6.4.1 Tensão elástica
A seção 3.3 do capítulo 3 descreve o conceito de tensão. Mas a aplicação prática em uma
simulação consistente de elasticidade requer a construção de um tensor não linear, que contem-
ple não só a deformação linear, mas também a angular. A técnica de Müller et al. (2004) utiliza
o tensor Green-Saint-Venant para o cálculo da tensão elástica, definida como:
ε = ∇u+∇uT +∇u∇uT (6.3)
51
Ou representando pela jacobiana J = I+∇uT :
ε = JJT − I (6.4)
Resultando em uma matriz 3×3 que descreve a tensão quadrática na abrangência do suporte
do phyxel em análise, sendo posteriormente utilizada para calcular o estresse de acordo com as
propriedades do material.
6.4.2 Tensão plástica
O método simula plasticidade como a absorção de parte da tensão elástica. Cada phyxel
deve guardar um tensor plástico e na etapa de cálculo do estresse operar sobre a diferença dos
tensores, com a tensão resultante definida por:
εresultante = εelastica− εplastica
Essa modelagem de plásticidade constroi-se sobre três parâmetros: max, yield e creep.
(GROSS; PFISTER, 2007)
Max representa o limite escalar para a quantidade de tensão que pode ser absorvida em forma
de deformação plástica. Ou seja, um phyxel só pode armazenar tensão plástica até um
valor máximo. Grandes valores para max causam o derretimento dos objetos sobre ten-
sões fortes o bastante. Um valor nulo de max significa que o objeto não possui comporta-
mento plástico.
Yield é o valor que define a tolerância a partir da qual a tensão elástica extravaza e é absorvida
como tensão plástica. Simula o limite de tensão que a estrutura molecular do corpo su-
porta até se romper e impedir que a elasticidade restaure a forma original. Se sofrer uma
tensão menor que o valor yield, nenhuma energia plástica é absorvida, se a tensão for
maior, então parte da deformação será absorvida e conservada no tensor plástico.
Creep é a taxa de tensão elástica que é absorvida e transformada em tensão plástica. Quando
a tensão elástica supera o valor max, soma-se à tensão plástica o valor da tensão elástica
multiplicado pelo creep.
Na simulação efetua-se os seguintes passos para calcular a absorção elástica e tensão resultante:
1. εresultante = εelastica− εplastica
52
2. Se ‖ εresultante ‖> yield
• εplastico = εplastico + εelastico · creep
3. Se ‖ εplastica ‖> max
• εplastico = εplastico · max‖εplastica‖
Assim atualizando a cada iteração a tensão plástica absorvida e calculando o estresse com base
na tensão resultante.
6.5 Calculando o estresse
Após o cálculo da tensão elástica resultante, a obtenção do tensor de estresse é trivial.
Seguindo as definições vistas na seção 3.4 do capítulo 3, e com o objetivo de simular materiais
isotrópicos, o cálculo de estresse constitui diretamente em:
σ = Cε
A matriz C é construida com base em dois parâmetros que descrevem o comportamento
elástico do material: módulo Young e a taxa de Poisson.
Young é a propriedade do corpo que mapeia linearmente o estresse com a força elástica resul-
tante.
Poisson é a taxa de preservação de volume.
Os dois valores são obtidos empiricamente sobre cada material.
Tendo os valores, constroi-se a matriz C (GROSS; PFISTER, 2007) como:
E(1+ v)(1−2v)
1− v v v 0 0 0
v 1− v v 0 0 0
v v 1− v 0 0 0
0 0 0 1−2v 0 0
0 0 0 0 1−2v 0
0 0 0 0 0 1−2v
(6.5)
53
Por C ser uma matriz 6× 6, a multiplicação não é feita pela matriz de tensão, e sim pelo
vetor coluna dos termos necessários, devido à já mencionada simetria de ε em relação à diagonal
principal. O vetor com as componentes do estresse fica definido como:
σ11
σ22
σ33
σ12
σ13
σ23
=
E(1+ v)(1−2v)
1− v v v 0 0 0
v 1− v v 0 0 0
v v 1− v 0 0 0
0 0 0 1−2v 0 0
0 0 0 0 1−2v 0
0 0 0 0 0 1−2v
ε11
ε22
ε33
ε12
ε13
ε23
(6.6)
A partir do tensor de estresse pode-se obter as forças elásticas resultantes.
6.6 Calculando as forças resultantes
O cálculo das forças resultantes da tensão e estresse se baseia na densidade da energia de
tensão acumulada em torno do phyxel, esta é definida por Müller et al. (2004) como:
U =12(ε ·σ)
E a força por unidade de volume como o gradiente negativo da energia de tensão em relação
ao deslocamento do phyxel. Resultando em:
f =−∇uU =−12
∇u(ε ·Cε) =−σ∇uε
Mas as forças resultantes são as forças aplicadas a cada um dos vizinhos, e o oposto da
soma destas atuando no phyxel sobre o qual está se iterando, devido ao princípio da ação e
reação. Então o cálculo efetivo de tais vetores de força é assim definido:
fi =−2viJiσidi = Fedi (6.7)
f j =−2viJiσid j = Fed j
Onde v é o volume do suporte do phyxel e:
54
di = A−1
(−∑
jxi jwi j
)(6.8)
d j = A−1xi jwi j
Tendo registrado a força recebida por todos os phyxels, pode-se prosseguir para a integração
temporal que irá gerar o movimento resultante da resposta elástica e plástica da simulação até
então.
6.7 Integração temporal
Ao final do cálculo de todas as forças elásticas cada phyxel possuirá um vetor de força
resultante acumulada da influência de todos os vizinhos sobre ele. Para simular o movimento
decorrente destas forças e prosseguir com a simulação, deve-se calcular a aceleração a gerada
pelas forças, a velocidade v gerada pela aceleração, e o deslocamento à uma nova posição
x gerado pela velocidade. Isto resume o conceito de integração explícita, especificamente a
integração numérica de Euler (SOUZA, 2011), que consiste na aplicação das leis de mecânica
clássica sobre um determinado intervalo de tempo:
a =ft+4t
mi
vt+4t = vt +a4t
xt+4t = xt +vt4t
Sua aplicação é simples e de baixo custo computacional, mas a simulação se torna bastante
sensível ao intervalo de tempo. Um 4t grande significa poucas iterações por segundo e em
consequência uma simulação mais rápida. Contudo, a simulação pode se tornar altamente in-
stável. Com grandes saltos no tempo, forças que deveriam atuar apenas em pequenos intervalos
acabam levando as partículas a se afastarem causando grandes deslocamentos. Isso resultaria
em forças ainda maiores, que no próximo passo de integração intensificariam o processo até
que o objeto perdesse coesão entre os phyxels, ou estourando a capacidade da notação ponto
flutuante da máquina.
Um integrador temporal explícito que ameniza este problema é o leap-frog (NOBREGA,
2010), ainda mantendo-se computacionalmente simples. Sua ideia é fragmentar a atualização
55
das informações de movimento da partícula, intercalando a atualização da velocidade e da acel-
eração, utilizando subintervalos na metade do4t para calcular a velocidade. Por isso seu nome
conotando o movimento em saltos de um sapo, pois a aceleração e velocidade atualizam-se em
saltos distintos no tempo, que acontecem da seguinte forma:
a =ft+4t
mi
vt+4t2= vt−4t
2+a4t (6.9)
xt+4t = xt +vt−4t24t
Que pode ser unificado e realizado apenas uma vez a cada intervalo de integração, aumen-
tando a estabilidade da integração sem elevar muito o custo computacional (HUT; MAKINO;
MCMILLAN, 1995), onde:
a =ft
mi
vt+4t = vt +12(at+4t +at
)4t (6.10)
xt+4t = xt +vt4t +12
at (4t)2
Os autores informam ainda sobre a possibilidade de integração implícita que considera o
estado de todas as partículas, e que exige a resolução de um sistema linear. Eles fornecem, in-
clusive, as matrizes necessárias já deduzidas que facilitam o processo que garantiria estabilidade
incondicional, mas com elevado custo computacional (MÜLLER et al., 2004).
6.8 Visualização
Assim como o interior é representado por partículas, nesta técnica a superfície é amostrada
por pontos chamados de surfels. Semelhantes aos phyxels, eles são partículas com coordenadas
de material, e durante o movimento tem posições globais que definem seu deslocamento. Pos-
suem uma vizinhança e uma função de peso em seu suporte. Diferem, porém, nas características
físicas que são negligenciadas. A massa, energia, força, aceleração e velocidade não entram em
consideração. Os surfels possuem apenas informações da vizinhança com os phyxels, sua nor-
mal em relação à superfície e sua posição.
56
A técnica baseia-se fortemente em um tipo de renderização diferente da convencional malha
de polígonos. Os surfels são movidos de acordo com o andamento da simulação, e são unifi-
cados e interpolados a partir da sua posição e da direção da sua normal, gerando a superfície
de aparência contínua no final da renderização (GROSS; PFISTER, 2007). Com ele, surge a
possibilidade de inserir dinamismo realista nos casos de grandes deformações, em especial do
derretimento. Quando um material é deformado de maneira que partes se expandam grande-
mente, os surfels se tornam esparsos e a superfície perde detalhamento. Então quando forem
expandidas ou quando novas áreas são criadas pela separação de partes do corpo, deve-se criar
novos surfels para preencher estes espaços. O inverso acontece quando o derretimento plástico
mergir duas faces da superfície, caso em que os surfels destas faces ficariam internos e poderiam
ser descartados.
Imparcial ao método de renderizar, a informação contida em cada surfel continua sendo
chave na representação visual deste método. Esta deve, então, ser atualizada de maneira coer-
ente de acordo com o movimento interno dos phyxels. Müller et al. (2004) expõe três formas
de realizar isso: uma pela consideração direta do deslocamento dos phyxels, uma pelo cálculo
da superfície implícita e uma abordagem mista. Este trabalho apresenta e utiliza apenas a mais
simples e direta delas, a atualização pelo deslocamento, explicada a seguir.
6.8.1 Atualização dos surfels pelo deslocamento
A ideia nesta abordagem é fazer com que os surfels acompanhem o movimento dos phyxels.
Onde o deslocamento de um surfel é o resultado combinado do deslocamento dos seus phyxels
vizinhos. Isto precisa ser feito com base em um campo vetorial de deslocamento suave que
defina a combinação do deslocamento de cada phyxel, ou seja, precisa-se das derivadas que for-
mam o gradiente do deslocamento ∇u. Isto já está disponível pela aplicação do MLS realizada
anteriormente e pode ser reutilizado na etapa de atualização da visualização.
Com base no conceito de partição de domínio que unifica a influência de todos os suportes
usando uma função de peso w em relação a cada distância di, o deslocamento ui de cada phyxel
da vizinhança e seu gradiente ∇ui, o deslocamento us f l de cada surfel pode ser calculado como:
us f l =1
∑i w(di)∑
iw(di)
(ui +∇uT
i(xs f l−xi
))(6.11)
Que é o deslocamento médio de cada phyxel vizinho projetado sobre o gradiente do desloca-
mento em relação às distâncias originais, ponderado pelas influências de seus pesos calculados
pela função Gaussiana:
57
w(di) = w(‖ xs f l−xi ‖
)= e
−d2
r2
Aplicando isto após a integração temporal dos phyxels estabelecemos uma atualização con-
sistente da superfície reposicionando cada surfel de acordo com seus vizinhos.
58
7 Implementação
Como forma de avaliação direta da eficiência e aplicabilidade da técnica de simulação elás-
tica e plástica baseada em pontos (MÜLLER et al., 2004), realizou-se a implementação com-
putacional do método. O fruto disso é um protótipo interativo que carrega objetos tridimen-
sionais, gera e configura as partículas físicas para a simulação, inicia e controla a iteração e
renderiza o resultado como uma malha de polígonos.
Nele o objeto em simulação está contido em um paralelepípedo onde há um objeto sólido
estático. O corpo deformável sofre colisão tanto com a caixa que envolve o cenário e com o
objeto estático. Há ainda a interferência por clique e arrasto do mouse sobre o objeto, o que
aplica uma força externa sobre o mesmo. Adicionalmente, há entradas para inserção e edição
dos valores que representam os atributos da simulação.
As seções a seguir descrevem os detalhes e estrutura de sua implementação, bem como
considerações práticas sobre o uso do método.
7.1 Recursos e bibliotecas
O protótipo resultante desta pesquisa foi feito na linguagem C++, com uso de recursos
de renderização, simulação física, técnicas e abstrações matemáticas disponíveis na biblioteca
C3DE, Cyclops 3D Enviroment (SILVA et al., 2009). Ela foi desenvolvida para servir como
arcabouço de métodos para renderização tridimensional na visualização de dados médicos e
como framework de suporte para simulações físicas, ainda fornecendo interface para periféricos
de interação humana como aparelhos dotados de acelerômetro. Sua arquitetura apresentando
seus componentes é vista na figura 7.1.
As classes desenvolvidas na implementação deste trabalho foram agregadas a essa bib-
lioteca, que passa a dispor dos algoritmos pesquisados para a técnica baseada em pontos. Per-
mitindo uso em futuras aplicações de simulação física de corpos deformáveis.
59
Figura 7.1: Visão superficial dos componentes que compõe o C3DE. (SILVA et al., 2009)
7.2 Arquitetura do protótipo
A biblioteca C3DE contempla a simulação de fluidos por SPH, o que já fornece abstração
de partícula bastante próxima da necessária para a construção da classe Phyxel, com atributos
de movimento: força, aceleração, velocidade, posição.
Outro fator importante é a utilização da workstation desenvolvida sobre esta biblioteca, que
atualmente suporta o acoplamento de plugins. Portanto, o protótipo foi desenvolvido como um
plugin que faz uso dos recursos de janela, renderização e interatividade da visualização.
Sua interface gráfica com o usuário gera eventos e fornece as variáveis da simulação. Uma
classe controla a comunicação entre a GUI e a classe responsável pela coordenação das etapas
da simulação, o WorldSimulator. Este é responsável por dar sequência à chamada de operações
de acordo com o definido no método, bem como a integração temporal realizada ao final de
cada iteração e receber as ações do usuário transmitidas pelo controlador na forma de forças
internas ou alteração de parâmetros.
Esta arquitetura foi planejada para ser simples, com o intuito de mostrar que pode-se re-
alizar simulações com o método diretamente sobre uma aplicação básica, desde que estejam
implementados os cálculos necessários em cada etapa da técnica. Ou seja, concentrou-se o
problema no algoritmo em si.
7.3 Carregamento de objetos
A principal motivação de realizar pesquisa em simulação de corpos deformáveis sem malha
é livrar-se da complexidade de representar os objetos computacionalmente. Com interesse
60
nisso, a forma de carregar objetos no simulador também partiu do princípio da simplicidade.
Pela necessidade do método, deve-se ter um volume que represente um corpo, para isto
usou-se modelos descritos em arquivos OBJ, compostos por malhas de triângulos.
Em seguida, este corpo deve ser preenchido com uma nuvem de phyxels. O método tolera
qualquer distribuição das partículas, mas para maior controle e observação dos efeitos esta
implementação usou uma forma regular para disposição dos phyxels. Isto foi feito criando
um algoritmo simples que vasculha a bounding box do objeto marchando segmentos de reta e
testando sua colisão contra a malha de triângulos do modelo usando o teste de intersecção por
árvore AABB (NOBREGA; B.; WANGENHEIM, 2008) previamente disponível no C3DE. O
algoritmo 1 descreve o processo.
Algoritmo 1 Algoritmo de preenchimento de phyxels no interior de um modelo.
Dado a c a i x a C , segmento de r e t a R e número de c o l i s õ e s N, f a ç a :P o s i c i o n e R num dos v é r t i c e s de C .Enquanto R não e s t i v e r além do topo da c a i x a f a ç a :
P o s i c i o n e R na d i r e i t a da c a i x a .Enquanto R não e s t i v e r além da l a t e r a l e s q u e r d a da c a i x a f a ç a :
N := 0 .P o s i c i o n e R na f r e n t e da c a i x a .Enquanto R não e s t i v e r além do fundo da c a i x a f a ç a :
Some o número de c o l i s õ e s com a malha em N.Se N f o r impar , c r i e um Phyxe l no c e n t r o de R .Mova R um p a s s o ao fundo .
Mova R um p a s s o à e s q u e r d a .Mova R um p a s s o acima .
Com ele qualquer malha de triângulos pode ser preenchida desde que ela seja completa-
mente fechada, sem buracos. O resultado, como exemplificado pela figura 7.2, não é preciso, e
depende do tamanho do segmento de reta utilizado. Passos menores resultam em maior número
de partículas distribuidas e melhor a aproximação da forma do objeto. Outro problema que
surgiu durante os testes foi o de falhas na contagem de colisões quando o objeto tem triângulos
que pertencem às faces da caixa envoltória. Isto faz com que phyxels sejam colocados fora do
objeto até reiniciar a contagem de colisões.
61
0 0 1 1 1 2 2
Figura 7.2: Demonstração de resultado do algoritmo de preenchimento com phyxels. A primeira
linha, em destaque, exibe a contagem de colisões desta passada através da figura, com os phyxels
colocados nos pontos ondem o número é impar.
Quando o tamanho do passo for pequeno o bastante e o alinhamento das faces não causar
problemas na detecção de colisão, o corpo obterá a lista de phyxels que foram distribuídos
adequadamente dentro dele, e estes poderão ser visualizados no protótipo.
Figura 7.3: Do modelo em malha de triângulos visto na esquerda, obteu-se a nuvem de partícu-
las exibida na direita em resultado do algoritmo de preenchimento.
7.4 Simulação
Uma vez carregado algum objeto e gerada sua nuvem de phyxels, pode-se iniciar os cálculos
preliminares da simulação. Cada phyxel deve ter sua vizinhança definida e ter seus parâmetros
configurados da maneira definida pelo método.
62
Além da definição da vizinhança, algumas informações podem ser calculadas e reutilizadas
ao longo de toda a simulação. Em cada phyxel, para cada vizinho, o método utiliza a cada iter-
ação o vetor de distância entre as posições iniciais, ou coordenadas de material. Para otimizar o
processo os objetos da classe Phyxel mantém uma lista com as distâncias calculadas assim que
cada vizinho é adicionado. Adicionalmente, sobre essas distâncias já se aplica a função de peso
para evitar a repetição desta operação.
O quadrado ponderado da distância original a cada vizinho é somado ao atributo que contém
a matriz A, definida em 6.1. Quando todos os vizinhos forem adicionados, pode-se realizar a
inversão de A com uso do SVD.
7.4.1 Configuração dos Phyxels
Müller et al. (2004) sugere a forma de configuração dos phyxels. Mas para maior com-
preensão do efeito destes valores de raio, massa e do número de vizinhos de cada partícula, a
implementação deste trabalho configura o raio do phyxel como a maior distância entre seus N
vizinhos mais próximos. Ressaltando que estes são obtidos por força bruta, já que isto só acon-
tece uma vez ao longo de todo o processo. Em aplicações que vão além do escopo de teste e
avaliação deste trabalho, sugere-se utilizar hashing espacial para encontrar vizinhos mais próx-
imos (TESCHNER et al., 2003). Desta forma o número de vizinhos oscila próximo de N, o que
permite a observação dos efeitos da variação do número de partículas conectadas ao suporte.
Observou-se que um número pequeno de vizinhos gera grande instabilidade na simulação,
provavelmente por somar os efeitos de conferir massa muito pequena e ter a matriz A mal condi-
cionada. Em alguns objetos, especialmente os com reentrâncias, um número muito grande de
vizinhos causou interferências globais indesejadas, pois seções delgadas e distantes acabavam
interferindo umas nas outras.
7.4.2 Aplicando SVD no cálculo da inversa de A
A classe Phyxel possui um método que realiza a inversão matricial por SVD da forma de-
scrita na subseção 6.3.2. Ele recebe o parâmetro que decide se o auto-valor é pequeno demais.
Após alguns testes verificou-se que para evitar a formação de componentes de força muito
grandes em casos de mal condicionamento, o que causaria instabilidade com a integração ex-
plícita, bons valores para ε encontram-se em torno de 1×10−3, de acordo com a observação do
estresse gerado em diversas deformações testadas. Valores menores que esse geram forças de
quarta ordem de magnitude, causando deslocamentos anômalos ao aplicar a integração tempo-
63
ral, e com isso um efeito de reforço que aumenta as forças resultantes até que o corpo perde a
coerência entre as partículas.
7.4.3 Colisão com superfícies estáticas
Para observar a dinâmica da deformação de um corpo é ideal que ocorram colisões que
causem grandes deslocamentos. Utilizou-se novamente os recursos de colisão por árvore AABB
já mencionados, desta vez para causar a colisão com a caixa que serve como delimitador do
mundo simulado. Além dela, inseriu-se um modelo estático, também com a mesma estrutura
de árvore, contra o qual pode-se chocar o objeto e observar a deformação resultante.
Nessas colisões, os phyxels são chocados contra as faces dos objetos estáticos e tem seu
movimento refletido de acordo com a normal do contato. Idealmente implementaria-se esta
colisão com os surfels para então distribuir aos phyxels, mas por razões de simplicidade e veri-
ficação mais imediata da dinâmica da deformação optou-se pelo simples contato com os phyxels
e apenas carregar os surfels junto do movimento.
Figura 7.4: Resultado da simulação de colisão entre modelo deformável e uma esfera estática.
A consequência da colisão é uma grande deformação no objeto, cujas forças elásticas atuam
tentando reestabelecer a forma original do corpo. O resultado deste efeito pode ser visto na
figura 7.4 que exibe dois quadros da simulação em tempo real ocorrendo em um caso de colisão
com um corpo estático.
64
7.4.4 Aplicação de forças externas
O protótipo permite que usuário aplique força externa sobre o corpo deformável usando o
clique e arrasto do mouse.
Quando o botão é pressionado, projeta-se uma linha saindo da posição do mouse em direção
ao fundo da cena. Verifica-se qual surfel essa linha atinge primeiro, e então seleciona-o.
Tendo um surfel selecionado, o movimento do mouse é rastreado enquanto o botão estiver
pressionado. Um vetor de força é criado de acordo com a direção e prolongamento do movi-
mento feito pelo usuário. Este vetor é, então, aplicado à todos os phyxels vizinhos do surfel
selecionado. Isto reproduz a ideia de uma pinça, podendo puxar e empurrar o corpo para obser-
var as consequências, como a figura 7.5 demonstra.
Figura 7.5: Três quadros que exibem o resultado da ferramenta de pinçar, a caixa vermelha
marca o local do seletor.
7.4.5 Atualização da visualização
Diferente do trabalho apresentado por Müller et al. (2004), esta implementação faz uso da
renderização convencional de objetos compostos por faces poligonais. Para preservar o conceito
do método, fez-se correlação entre os vértices do modelo com os surfels do objeto.
Quando um modelo é carregado cria-se um surfel para cada vértice, que depois tem suas
vizinhanças configuradas da mesma forma que os phyxels. Com exceção do fato que um surfel
tem phyxels como vizinhos, e não outros surfels, pois eles não interferem entre si.
No decorrer da simulação eles tem sua posição atualizada de acordo com o deslocamento
65
dos phyxels. Isso modifica as faces do modelo para o novo estado da animação criando a ideia
de movimento e exibindo diretamente o resultado da simulação.
Mas não basta atualizar a informação dos vértices lendo a posição dos surfels. A iluminação
aplicada na etapa de renderização depende da informação das normais dos vértices para gerar
efeito convincente de profundidade. Sendo assim, ao final da atualização da posição dos vértices
deve-se recalcular e atualizar suas normais. Isto é feito com produto cruzado com todos os pares
de arrestas conectadas pelo vértice sobre o qual está se iterando para atualizar sua normal. O
produto cruzado fornece um vetor perpendicular a cada uma dessas arrestas, e a soma desses
vetores gera a direção média, que é então normalizada e utilizada para a etapa de iluminação.
Como recurso extra na visualização, pode-se optar entre exibir o objeto opaco, ou com
transparência parcial. Na segunda opção também são renderizadas na cena esferas que marcam
a posição de cada phyxel ao longo da simulação. Isto permite observar seu comportamento
e como ocorre a dinâmica elástica e plástica, incluindo o resultado do algoritmo de preenchi-
mento, como já demonstrado na figura 7.3.
7.5 Resultados
Com o protótipo concluído, testes visuais e de desempenho foram realizados para analisar o
potencial da técnica. Diversos modelos em malha de triângulos foram carregados pela aplicação
e simulados com variações de parâmetros físicos e de distribuição das partículas.
Para valores intermediários de módulo Young, no geral entre 1K e 40K, o comportamento
elástico seguiu o esperado, gerando animações coerentes das deformações sofridas pelos objetos
ao colidirem com superfícies estáticas e pelo resultado da aplicação de força com a ferramenta
de pinçar e arrastar.
A simulação se mostrou bastante sensível à variação da quantidade de phyxels preenchendo
o interior dos corpos. Aproximações rudes puderam ser feitas para formas simples usando
poucas partículas, o que resultou numa simulação rápida que poderia ser usada em aplicações
que descartem grande realismo. Aumentando a quantidade de phyxels, formas mais complexas
puderam ser aproximadas, como rostos e orgãos humanos. Aqui o desempenho cai consider-
avelmente, mas o resultado visual seguiu o esperado pelo princípio da partição de domínio,
refinando o comportamento elástico do material.
A variação do número de partículas causou, ainda, outro efeito colateral: a quantidade de
massa distribuida em cada porção modifica conforme a granularidade muda. Portanto, deve-se
66
Phyxels Surfels Forças Atualização56 482 0.20956014 ms 0.1355661 ms
106 7136 0.33372935 ms 2.1721624 ms498 482 1.7688644 ms 0.15197148 ms850 7136 3.0448024 ms 2.35544 ms1760 482 6.8089768 ms 0.23258457 ms2844 7136 9.4171316 ms 2.486963 ms
Tabela 7.1: Resultados em milissegundos da simulação de uma esfera com 482 surfels e domodelo da cabeça de Stalin com 7136 surfels, demonstrando a relação linear entre o número dephyxels e o cálculo das forças elásticas efetuado pelo protótipo rodando em uma máquina comprocessador Intel de quatro núcleos de 3,20GHz de frequência e 8 MB de memória RAM.
usar um fator de escala que equilibre a densidade do material para que não seja necessário alterar
o módulo Young para atingir o mesmo comportamento de elasticidade. Ou seja, menor número
de partículas precisaria de um menor valor no módulo Young para simular um material mais
rígido do que o mesmo corpo com mais partículas. Isso deve ser feito na etapa de configuração
dos pesos de cada phyxel, onde sua massa é calculada.
A variação de partículas afeta diretamente a velocidade da simulação, já que o algoritmo
depende linearmente da quantidade de phyxels e do número de vizinhos que cada um possui,
correspondendo à complexidade O(N +NM), onde N é o número de partículas e M o número
de vizinhos, como pode ser verificado na tabela 7.1.
A possibilidade de variar a quantidade de vizinhos demonstrou que formas mais complexas
e com grande número de partículas em seu interior requererem maior número de vizinhos por
phyxel, em parte pelo condicionamento da matriz de momento, e também pela distribuição da
massa que depende desse valor. Para formas mais simples, um número menor, com um mínimo
de 7 vizinhos, pode ser usado para leve aumento na velocidade da simulação.
Quanto à integração explícita, usou-se intervalo de tempo de 1× 10−3, valor que permite
a configuração de materiais com módulo Young de até 40K de maneira estável para a maioria
dos casos, com uma quantidade mediana de partículas. Em corpos gerados com mais de 500
phyxels, a soma do tempo de cálculo das forças e atualização da malha aproxima-se de 10 milis-
segundos, indo além da capacidade de realizar um quadro de simulação dentro do intervalo de
tempo da integração, fazendo com que o retorno visual seja mais lento do que deveria. Com a
adição de técnicas de detecção de colisão entre corpos deformáveis ou outras funcionalidades,
a soma dos tempos elevaria-se ainda mais. Deve-se, então, explorar o amplo potencial de par-
alelismo da técnica para alcançar desempenho satisfatório para aplicações interativas, já que
cada partícula pode ter as forças resultantes calculadas individualmente.
67
8 Conclusão
O grande contraste da classe de métodos livres de malha e os dependentes de malha se
mostrou evidente pela trivialidade como corpos de formas complexas puderam ser representa-
dos por nuvens de pontos criadas com algoritmo simples sem qualquer preocupação estrutural.
Os casos estáveis tiveram comportamento realista sob análise visual, sem apresentar artefatos
anômalos de deformação angular, mostrando que a aplicação do tensor quadrático oferece re-
spostas elásticas consistentes, mesmo sob grandes deformações e variações da forma original.
Contudo, a frequente ocorrência casos instáveis deixou evidente que para uma aplicação efetiva
do método a distribuição dos phyxels não pode ser descuidada. Corpos complexos, com partes
delgadas ou concavidades, deveriam ter suas distintas seções analisadas individualmente. Para
isso elaboraria-se um algoritmo de preenchimento mais específico, criando distribuições que
favoreçam a estabilidade numérica na integração explícita. Portanto, a primeira indicação é de
que a independência de malha não exime o método de um refinamento extra na construção das
representações dos objetos.
Quanto ao desempenho alcançado, sua variabilidade é grande , assim como a dependên-
cia da estabilidade. A implementação de um integrador explícito buscou realizar a troca de
estabilidade por eficiência. Mas uma aplicação interativa não serve se for rápida, mas perder
a coerência e desmantelar a simulação. Então, para uma aplicação mais robusta precisaria-se
reduzir o passo de integração, ou utilizar integrador implícito sugerido por Müller et al. (2004).
Mas ambas afetam o desempenho, um passo de integração menor significa mais atualizações
físicas por segundo e em consequência mais instruções realizadas. Outro ponto delicado é da
quantidade de partículas utilizadas, ligadas diretamente ao abordado sobre a discretização do
contínuo. Resultados mais realistas requerem maior número de phyxels, resultando diretamente
em maior número de operações a cada iteração. Então a aplicação deste método em protótipo
interativo deve ponderar sobre o nível de realismo visual desejado.
O grande benefício que pode ser explorado é o potencial de se usar paralelismo nos cálculos
efetuados em cada iteração. Pois cada phyxel pode ter as forças elásticas calculadas independen-
temente dos outros. Sob ligeira análise informal, o único ponto crítico é o acûmulo das forças,
68
já que resultados indepentendes calculados em phyxels distintintos conflitariam ao agregar as
forças resultantes sobre vizinhos em comum. Depois de calcular os resultados e controlar a con-
corrência no registro das forças, a integração explícita poderia acontecer também em paralelo
para cada phyxel. O mesmo se aplica à atualização dos surfels. Então, no quesito eficiência
e estabilidade, recomenda-se futuro de uso de hardware realize as operações em paralelo, em
especial utilizando GPUs (FARIAS et al., 2008).
Quanto à validação numérica, estudos sobre a consistência da simulação em aplicações de
análise científica e técnica, como cálculo estrutural para predição de comportamento, devem
ser realizados antes que se possa afirmar algo sobre o método. Apesar de ser derivado de
formulações consolidadas da mecânica contínua, seus resultados são, a princípio, adequados
apenas para aplicações visuais.
A renderização utilizada no protótipo foi uma adaptação distinta da visualização proposta
pelo método original. Devido a ela, grandes deformações e derretimento do corpo tornam
incoerente o modelo poligonal da superfície, com faces invertidas e entrelaçadas desvirtuando
a visualização do objeto. Portanto, para melhor aproveitamento dos resultados da simulação
física recomenda-se implementar a renderização utilizada por Müller et al. (2004) que pode
lidar com grandes deformações. Assim descartando efetivamente o uso de malha no processo.
Mesmo com a forma alternativa de renderização, os resultados das simulações realizadas
no protótipo fruto da aplicação do método baseado em pontos exibiram grande apelo visual.
O desenvolvimento da função de interação com efeito de pinça, realizada também de maneira
simples e direta, demonstra o potencial do método para animações interativas. Unindo este
recurso à colisão com objetos estáticos, abre-se a possibilidade de aplicações na área médica,
como exemplo: avaliação da posição de próteses rígidas no interior de estruturas anatômicas
com comportamento elástico, como músculos e pele, e utilizar o efeito visual da plásticidade
para criação de moldes virtuais.
Isto mostra que o uso de métodos sem malha é viável e oferece grandes possibilidades,
especialmente por ser uma linha de pesquisa em pleno desenvolvimento. Estudos futuros de-
vem consolidar ainda mais os benefícios oferecidos por esta classe, talvez permitindo que eles
suplantem os métodos baseados em malha.
69
Referências Bibliográficas
BELYTSCHKO, T.; LU, Y. Y.; GU, L. Element-free galerkin methods. International Journalfor Numerical Methods in Engineering, John Wiley and Sons, Ltd, v. 37, n. 2, p. 229–256,1994. ISSN 1097-0207. Disponível em: <http://dx.doi.org/10.1002/nme.1620370205>.
BELYTSCHKO, T.; TABBARA, M. Dynamic fracture using element-free galerkinmethods. International Journal for Numerical Methods in Engineering, John Wileyand Sons, Ltd, v. 39, n. 6, p. 923–938, 1996. ISSN 1097-0207. Disponível em:<http://dx.doi.org/10.1002/(SICI)1097-0207(19960330)39:6<923::AID-NME887>3.0.CO;2-W>.
BOWEN, R. Introduction to Continuum Mechanics for Engineers: Revised Edition. DoverPublications, 2010. (Dover Civil and Mechanical Engineering Series). ISBN 9780486474601.Disponível em: <http://books.google.com.br/books?id=S8cWPwAACAAJ>.
CLKER. Disponível em: <http://www.clker.com/clipart-wireframe-3d-mesh-primitives.html>.
FARIAS, T. de; ALMEIDA, M.; TEIXEIRA, J.; TEICHRIEB, V.; KELNER, J. A highperformance massively parallel approach for real time deformable body physics simulation.In: Computer Architecture and High Performance Computing, 2008. SBAC-PAD ’08. 20thInternational Symposium on. [S.l.: s.n.], 2008. p. 45 –52. ISSN 1550-6533.
FELIPPA, C. A. Introduction to Finite Element Method. University of Colorado Boulder,Colorado 80309-0429, USA: University of Colorado, 2001.
FOTEINOS, P.; CHRISOCHOIDES, N. High-quality multi-tissue mesh generation for finiteelement analysis. In: MeshMed, Workshop on Mesh Processing in Medical Image Analysis(MICCAI). [S.l.: s.n.], 2011. p. 18–28.
GINGOLD, R. A.; MONAGHAN, J. J. Smoothed particle hydrodynamics - theory andapplication to non-spherical stars. Monthly Notices of the Royal Astronomical Society, vol. 181,p. 375-389., November 1977.
GROSS, M.; PFISTER, H. Point-Based Graphics. San Francisco, CA, USA: MorganKaufmann Publishers Inc., 2007. ISBN 0123706041, 9780080548821.
HEINBOCKEL, J. H. Introduction to Tensor Calculus and Continuum Mechanics. [s.n.], 2001.ISBN 978-1-55369-133-4. Disponível em: <http://www.math.odu.edu/ jhh/counter2.html>.
HUNTER, P.; PULLAN, A. FEM/BEM Notes. [S.l.]: University of Auckland New Zealand,2001.
HUT, P.; MAKINO, J.; MCMILLAN, S. Building a better leapfrog. Astrophysical Journal,Part 2 - Letters (ISSN 0004-637X), vol. 443, no. 2, p. L93-L96, 1995.
70
IDELSOHN, S. R.; ONATE, E. To mesh or not to mesh. that is the question. Com-puter Methods in Applied Mechanics and Engineering, v. 195, 2006. ISSN 0045-7825. <ce:title>John H. Argyris Memorial Issue. Part I</ce:title>. Disponível em:<http://www.sciencedirect.com/science/article/pii/S0045782505005013>.
KITCHENHAM, B. Procedures for performing systematic reviews. [S.l.], 2004.
LANCASTER, P.; SALKAUSKAS, K. Surfaces Generated by Moving Least Squares Methods.Mathematics of Computation, American Mathematical Society, v. 37, n. 155, p. 141–158,1981. ISSN 00255718. Disponível em: <http://dx.doi.org/10.2307/2007507>.
LEVIN, D. The approximation power of moving least-squares. Mathematics of Computation,v. 67, p. 1517–1531, 1998.
MOSEGAARD, J. Cardiac Surgery Simulation - Graphics Hardware meets CongenitalHeart Disease. Tese (Doutorado) — Department of Computer Science, University of Aarhus,Denmark, oct 2006.
MÜLLER, M.; HEIDELBERGER, B.; TESCHNER, M.; GROSS, M. Meshlessdeformations based on shape matching. ACM Trans. Graph., ACM, New York,NY, USA, v. 24, n. 3, p. 471–478, jul 2005. ISSN 0730-0301. Disponível em:<http://doi.acm.org/10.1145/1073204.1073216>.
MÜLLER, M.; KEISER, R.; NEALEN, A.; PAULY, M.; GROSS, M.; ALEXA, M. Pointbased animation of elastic, plastic and melting objects. In: Proceedings of the 2004 ACMSIGGRAPH/Eurographics symposium on Computer animation. Aire-la-Ville, Switzerland,Switzerland: Eurographics Association, 2004. (SCA ’04), p. 141–151. ISBN 3-905673-14-2.Disponível em: <http://dx.doi.org/10.1145/1028523.1028542>.
NAYROLES, B.; TOUZOT, G.; VILLON, P. Generalizing the finite element method:Diffuse approximation and diffuse elements. Computational Mechanics, Springer Berlin /Heidelberg, v. 10, p. 307–318, 1992. ISSN 0178-7675. 10.1007/BF00364252. Disponível em:<http://dx.doi.org/10.1007/BF00364252>.
NEALEN, A.; DARMSTADT, T. U. An as-short-as-possible introduction to the least squares,weighted least squares and moving least squares methods for scattered data approximation andinterpolation. URL: http://www. nealen. com/projects, 2004.
NEALEN, A.; MUELLER, M.; KEISER, R.; BOXERMAN, E.; CARLSON, M. Physicallybased deformable models in computer graphics. Blackwell Publishing, 2006. ISSN 0167-7055.Disponível em: <http://dx.doi.org/10.1111/j.1467-8659.2006.01000.x>.
NOBREGA, T.; B., D. D.; WANGENHEIM, A. von. Using metaprogramming func-tions to implement double-dispatch for collision. In: . [s.n.], 2008. Disponível em:<www.gpec.ucdb.br/sibgrapi2008/posters/47829.pdf>.
NOBREGA, T. d. H. C. Metodologia computacional para a simulação interativa de fluidosem estruturas tubulares. Dissertação (Mestrado) — Universidade Federal de Santa Catarina,2010.
PRESS, W. H.; TEUKOLSKY, S. A.; VETTERLING, W. T.; FLANNERY, B. P. Numericalrecipes in C (2nd ed.): the art of scientific computing. New York, NY, USA: CambridgeUniversity Press, 1992. ISBN 0-521-43108-5.
71
SILVA, A. F. B.; CARVALHO, D. D. B.; NOBREGA, T.; INáCIO, R. T.; WANGENHEIM,A. von. A multi-layered development framework for medical imaging applications. In:CBMS. IEEE, 2009. p. 1–4. ISBN 978-1-4244-4878-4. Disponível em: <http://dblp.uni-trier.de/db/conf/cbms/cbms2009.html>.
SMITH, C. On vertex-vertex systems and their use in geometric and biological modelling. Tese(Doutorado), Calgary, Alta., Canada, Canada, 2006. AAINR19574.
SOUZA, M. S. A rigid body physics engine for interactive applications. SBGames, 2011.
TESCHNER, M.; HEIDELBERGER, B.; MUELLER, M.; POMERANETS, D.; GROSS, M.Optimized spatial hashing for collision detection of deformable objects. In: . [S.l.: s.n.], 2003.p. 47–54.
ZACHOW, S.; GLADILINE, E.; HEGE, H. c.; DEUFLHARD, P. Finite-element simulationof soft tissue deformation. In: Computer Assisted Radiology and Surgery (CARS), ElsevierScience B.V. [S.l.: s.n.], 2000. p. 23–28.