Flow123d  3.9.1-60c7e5c
intersectionquadrature.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file intersection.hh
15  * @brief
16  */
17 
18 #ifndef INTERSECTION_HH_
19 #define INTERSECTION_HH_
20 
21 /*!
22  *
23  * Copyright (C) 2007 Technical University of Liberec. All rights reserved.
24  *
25  * Please make a following refer to Flow123d on your project site if you use the program for any purpose,
26  * especially for academic research:
27  * Flow123d, Research Centre: Advanced Remedial Technologies, Technical University of Liberec, Czech Republic
28  *
29  * This program is free software; you can redistribute it and/or modify it under the terms
30  * of the GNU General Public License version 3 as published by the Free Software Foundation.
31  *
32  * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
33  * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
34  * See the GNU General Public License for more details.
35  *
36  * You should have received a copy of the GNU General Public License along with this program; if not,
37  * write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 021110-1307, USA.
38  *
39  *
40  * $Id: neighbours.h 1055 2011-04-21 13:43:54Z jan.brezina $
41  * $Revision: 1055 $
42  * $LastChangedBy: jan.brezina $
43  * $LastChangedDate: 2011-04-21 15:43:54 +0200 (Thu, 21 Apr 2011) $
44  *
45  * @file
46  * @brief ???
47  *
48  */
49 
50 
51 #include <armadillo>
52 
53 #include "mesh/ngh/include/intersectionLocal.h"
54 
55 
56 
57 /**
58  * Navrh algoritmu pro hledani pruniku elementu dvou siti (libovlnych dimenzi)
59  * algoritmus postupuje od bodu pruniku pres usecky a polygony k mnohostenum
60  *
61  * Vstup: Sit1 dimenze d1 a Sit2 dimenze d2
62  * predpoladam d1<=d2
63  *
64  * 1) hladam body na hranici pruniku tj.
65  * Intersection<d> <d_e1,d_e2> .. prunik ma dimenzi d a pronikaji se simplexy dimenze d_e1, a d_e2
66  *
67  * Intersection<0><0,0> .. totozne vrcholy El<0>
68  * Intersection<0><0,1> a <1,0> .. vrchol jedne site lezi na hrane druhe site
69  * Intersection<0><0,n> a <n,0> .. vrchol lezi na El<n> druhe site
70  *
71  * Intersection<0><1,1> .. bodovy prusecik dvou usecek v rovine
72  * Intersection<0><1,2> a <2,1> ... prusecik hrany a trojuhelnika
73  * ... dalsi zvlastni pripady vcetne <0><3,3> .. tetrahedrony s vrcholem na povrchu druheho
74  *
75  * 2) liniove pruniky Intersection<1>:
76  * Intersection<1><1,1> .. usecky na spolecne primce
77  * Intersection<1><1,2> a <2,1>.. usecka v rovine trojuhelnika
78  * Intersection<1><1,3> a <3,1> .. usecka a tetrahedron
79  * Intersection<1><2,2> .. prusecik dvou trojuhelniku
80  * Intersection<1><2,3> a <3,2> .. trojuhelnik a hrana tetrahedronu
81  * ..
82  *
83  * ... doprcic je to fakt hodne moznosti a je otazka, zda je nutne je vsechny rozlisovat
84  *
85  * Algoritmus by mel probuhat takto:
86  *
87  * 1) Najdu vrchol V site 1 a element E site 2 aby V byl v E
88  * (to neni tak trivialni, pokud site nepokryvaji stejnou oblast ale snad by to slo hledat v
89  * pruniku obalovych boxu)
90  * 2) najdu pruseciky P_i hran z vrcholu V s povrchem E,
91  * konstruuju vsechny potrebne pruniky elementu majici vrchol V s elementem E
92  *
93  * Sousedni elementy spolu s hranami ktere do nich vedou ulozim do prioritni fronty.
94  *
95  * 3) Vyberu z prioritni fronty novy E, pricemz vyuzivam spositane pruseciky psislusne steny a okoli vrcholu V
96  * tj. jdu po hranach po kterych jsem do noveho lementu prisel a najdu vsechny hranove pruniky, pak konstuuju slozitejsi pruniky
97  * az mam vsechny pruniky s novym elementem ...
98  *
99  * ...
100  *
101  * Prioritni fronta by preferovala elementy do kterych jsem se nejvicekrat dostal, tim se snazim minimalizovat povrch projite oblasti.
102  * Je ale mozne, ze to algoritmus naopak zpomali, pokud je prioritni fronta log(n).
103  *
104  * Zpracovani jednoho elementu tedy zahrnuje
105  * 1) trasovani hran:
106  * pro hranu H: testuju hledam prusecik se ctyrstenem:
107  * ANO -> pamatuju si hranovy prunik a ke stene (resp. sousednimu elementu) kde hrana vychazi pridam vychozi hranu
108  * NE -> konci ve vrcholu, dalsi hrany vychazejici z vrcholu pridam na seznam hran vchazejicich do elementu
109  *
110  * 2) po nalezeni pruniku vsech hran, hledam pruniky vsech vchazejicich ploch:
111  * jedna plocha ma se vstupni stenou useckovy prunik na jehoz konci jsou:
112  * *vstupni hrana
113  * * okraj steny
114  * kazdopadne trasuju okraj plosneho pruniku pres povrch elementu, nebo po vstupnich hranach,
115  * pokud prunikova plocha obsahuje vrchol tvorim nove vstupni plochy ...
116  *
117  * 3) Podobne trasuju vchazejici objemy
118  *
119  * ?? lze nejak vyuzit pokud ma element vice vstupnich sten
120  * minimalne se da kontrolovat ...
121  *
122  *
123  * Struktura systemu pruniku do budoucna:
124  * 1) trida IntersectionManager, ma matici vektoru. Na poli A(i,j) je vektor lokalnich souradnic na elementu dimenze i (chodi od 1 do 3)
125  * pruniku dimenze j (chodi od 0 do 3 resp do 2 pokud nebudu chtit prekryvy siti stejne dimenze)
126  *
127  * 2) Jeden intersection objekt je pak iterator dvou elementu a dva indexy lokalnich souradnic v prislusnych vektorech.
128  *
129  * Prozatim to zjednodusime tak, ze vektory lokalnich souradnic budu alokovat zvlast a nebudu je zdruzovat
130  *
131  *
132  * Nakonec potrebuju pocitat integral pres prunik z nejake funkce f(phi_a(x), phi_b(x)), kde phi_a je bazova funkce na jednom elementu a phi_b na druhem.
133  * To budu delat numerickou kvadraturou, takze potrebuji zobrazit prunik na jednotkovy simplex. Pro uzel kvardatury x_i musim najit body a_i a b_i na
134  * referencnich elementech A a B. Tj potrebuju lokalni souradnice (to jsou souradnice na referencnich elementech) kvadraturnich bodu. V nic pak umim spocitat hodnotu bazovych funkci
135  * a pak i hodnotu funkce f.
136  *
137  * K tomu staci mit matici transformace pruniku na referencni element. Takze bych pro jednotlive dvojice element - prunik mel matici + posouvaci vektor z armadila.
138  *
139  *
140  *
141  */
142 
143 #include <iostream>
144 template <int spacedim> class ElementAccessor;
145 
147 public:
148  Intersection(ElementAccessor<3> ele_master, ElementAccessor<3> ele_slave,
149  const IntersectionLocal *isec);
150 
151  /// dimension of the master element
152  unsigned int master_dim();
153 
154  /// dimension of the slave element
155  unsigned int slave_dim();
156 
157  const Element * master_iter() const
158  {return const_cast<Element *>( (Element *)(master) );}
159  const Element * slave_iter() const
160  {return const_cast<Element *>( (Element *)(slave) );}
161 
162  arma::vec map_to_master(const arma::vec &point) const;
163  arma::vec map_to_slave(const arma::vec &point) const;
164  double intersection_true_size() const;
165 
166  ElementAccessor<3> master, slave; // master lower dimension
167 
168 private:
169  /// dimenze pruniku
170  unsigned int dim;
171 
172  /// matrix part of linear transform from reference element of intersection to reference element of master or slave
173  arma::Mat<double> master_map, slave_map;
174  /// shift vector of the linear transform
176  void intersection_point_to_vectors(const IntersectionPoint *point, arma::vec &vec1, arma::vec &vec2);
177 };
178 
179 
180 #endif /* INTERSECTION_HH_ */
Intersection::master_iter
const Element * master_iter() const
Definition: intersectionquadrature.hh:157
Intersection::slave
ElementAccessor< 3 > slave
Definition: intersectionquadrature.hh:166
Armor::vec
ArmaVec< double, N > vec
Definition: armor.hh:885
Intersection::master
ElementAccessor< 3 > master
Definition: intersectionquadrature.hh:166
ElementAccessor
Definition: dh_cell_accessor.hh:32
Intersection::slave_shift
arma::vec slave_shift
Definition: intersectionquadrature.hh:175
Intersection::master_map
arma::Mat< double > master_map
matrix part of linear transform from reference element of intersection to reference element of master...
Definition: intersectionquadrature.hh:173
Element
Definition: elements.h:40
Intersection::intersection_point_to_vectors
void intersection_point_to_vectors(const IntersectionPoint *point, arma::vec &vec1, arma::vec &vec2)
Definition: intersectionquadrature.cc:68
Intersection
Definition: intersectionquadrature.hh:146
Intersection::dim
unsigned int dim
dimenze pruniku
Definition: intersectionquadrature.hh:170
Intersection::intersection_true_size
double intersection_true_size() const
Definition: intersectionquadrature.cc:101
Intersection::slave_dim
unsigned int slave_dim()
dimension of the slave element
Definition: intersectionquadrature.cc:63
Intersection::master_shift
arma::vec master_shift
shift vector of the linear transform
Definition: intersectionquadrature.hh:175
Intersection::master_dim
unsigned int master_dim()
dimension of the master element
Definition: intersectionquadrature.cc:58
Intersection::map_to_master
arma::vec map_to_master(const arma::vec &point) const
Definition: intersectionquadrature.cc:80
IntersectionLocal
Class represents intersection of two elements.
Definition: inspect_elements_algorithm.hh:34
IntersectionPoint
Class represents an intersection point of simplex<N> and simplex<M>. It contains barycentric coordina...
Definition: intersection_local.hh:32
Intersection::map_to_slave
arma::vec map_to_slave(const arma::vec &point) const
Definition: intersectionquadrature.cc:91
Intersection::slave_map
arma::Mat< double > slave_map
Definition: intersectionquadrature.hh:173
Intersection::Intersection
Intersection(ElementAccessor< 3 > ele_master, ElementAccessor< 3 > ele_slave, const IntersectionLocal *isec)
Definition: intersectionquadrature.cc:28
Intersection::slave_iter
const Element * slave_iter() const
Definition: intersectionquadrature.hh:159