Flow123d  master-f44eb46
intersection_quadrature.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 "mesh/accessors.hh"
52 #include <armadillo>
53 //#include <iostream>
55 
56 //#include "mesh/ngh/include/intersectionLocal.h"
57 
58 
59 
60 /**
61  * Navrh algoritmu pro hledani pruniku elementu dvou siti (libovlnych dimenzi)
62  * algoritmus postupuje od bodu pruniku pres usecky a polygony k mnohostenum
63  *
64  * Vstup: Sit1 dimenze d1 a Sit2 dimenze d2
65  * predpoladam d1<=d2
66  *
67  * 1) hladam body na hranici pruniku tj.
68  * Intersection<d> <d_e1,d_e2> .. prunik ma dimenzi d a pronikaji se simplexy dimenze d_e1, a d_e2
69  *
70  * Intersection<0><0,0> .. totozne vrcholy El<0>
71  * Intersection<0><0,1> a <1,0> .. vrchol jedne site lezi na hrane druhe site
72  * Intersection<0><0,n> a <n,0> .. vrchol lezi na El<n> druhe site
73  *
74  * Intersection<0><1,1> .. bodovy prusecik dvou usecek v rovine
75  * Intersection<0><1,2> a <2,1> ... prusecik hrany a trojuhelnika
76  * ... dalsi zvlastni pripady vcetne <0><3,3> .. tetrahedrony s vrcholem na povrchu druheho
77  *
78  * 2) liniove pruniky Intersection<1>:
79  * Intersection<1><1,1> .. usecky na spolecne primce
80  * Intersection<1><1,2> a <2,1>.. usecka v rovine trojuhelnika
81  * Intersection<1><1,3> a <3,1> .. usecka a tetrahedron
82  * Intersection<1><2,2> .. prusecik dvou trojuhelniku
83  * Intersection<1><2,3> a <3,2> .. trojuhelnik a hrana tetrahedronu
84  * ..
85  *
86  * ... doprcic je to fakt hodne moznosti a je otazka, zda je nutne je vsechny rozlisovat
87  *
88  * Algoritmus by mel probuhat takto:
89  *
90  * 1) Najdu vrchol V site 1 a element E site 2 aby V byl v E
91  * (to neni tak trivialni, pokud site nepokryvaji stejnou oblast ale snad by to slo hledat v
92  * pruniku obalovych boxu)
93  * 2) najdu pruseciky P_i hran z vrcholu V s povrchem E,
94  * konstruuju vsechny potrebne pruniky elementu majici vrchol V s elementem E
95  *
96  * Sousedni elementy spolu s hranami ktere do nich vedou ulozim do prioritni fronty.
97  *
98  * 3) Vyberu z prioritni fronty novy E, pricemz vyuzivam spositane pruseciky psislusne steny a okoli vrcholu V
99  * tj. jdu po hranach po kterych jsem do noveho lementu prisel a najdu vsechny hranove pruniky, pak konstuuju slozitejsi pruniky
100  * az mam vsechny pruniky s novym elementem ...
101  *
102  * ...
103  *
104  * Prioritni fronta by preferovala elementy do kterych jsem se nejvicekrat dostal, tim se snazim minimalizovat povrch projite oblasti.
105  * Je ale mozne, ze to algoritmus naopak zpomali, pokud je prioritni fronta log(n).
106  *
107  * Zpracovani jednoho elementu tedy zahrnuje
108  * 1) trasovani hran:
109  * pro hranu H: testuju hledam prusecik se ctyrstenem:
110  * ANO -> pamatuju si hranovy prunik a ke stene (resp. sousednimu elementu) kde hrana vychazi pridam vychozi hranu
111  * NE -> konci ve vrcholu, dalsi hrany vychazejici z vrcholu pridam na seznam hran vchazejicich do elementu
112  *
113  * 2) po nalezeni pruniku vsech hran, hledam pruniky vsech vchazejicich ploch:
114  * jedna plocha ma se vstupni stenou useckovy prunik na jehoz konci jsou:
115  * *vstupni hrana
116  * * okraj steny
117  * kazdopadne trasuju okraj plosneho pruniku pres povrch elementu, nebo po vstupnich hranach,
118  * pokud prunikova plocha obsahuje vrchol tvorim nove vstupni plochy ...
119  *
120  * 3) Podobne trasuju vchazejici objemy
121  *
122  * ?? lze nejak vyuzit pokud ma element vice vstupnich sten
123  * minimalne se da kontrolovat ...
124  *
125  *
126  * Struktura systemu pruniku do budoucna:
127  * 1) trida IntersectionManager, ma matici vektoru. Na poli A(i,j) je vektor lokalnich souradnic na elementu dimenze i (chodi od 1 do 3)
128  * pruniku dimenze j (chodi od 0 do 3 resp do 2 pokud nebudu chtit prekryvy siti stejne dimenze)
129  *
130  * 2) Jeden intersection objekt je pak iterator dvou elementu a dva indexy lokalnich souradnic v prislusnych vektorech.
131  *
132  * Prozatim to zjednodusime tak, ze vektory lokalnich souradnic budu alokovat zvlast a nebudu je zdruzovat
133  *
134  *
135  * 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.
136  * 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
137  * 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
138  * a pak i hodnotu funkce f.
139  *
140  * K tomu staci mit matici transformace pruniku na referencni element. Takze bych pro jednotlive dvojice element - prunik mel matici + posouvaci vektor z armadila.
141  *
142  * TODO:
143  * - review previous notes
144  *
145  */
146 
147 
148 /**
149  * Quadrature class for mortar like P0 algorithm. Currently is considered to be constructed on fly from IntersectionLocal.
150  * If we find a way how to merge it with IntersectionQuadratureP1 we can save these new object instead IntersectionLocal
151  * or precompute these if it pays off.
152  */
154 public:
156  : mesh_(mesh)
157  {}
158 
159  /**
160  * Reinit quadrature measure. Returns true if quadrature is close to zero.
161  */
162  inline bool reinit(const IntersectionLocalBase *isec)
163  {
164  slave_idx_ = isec->bulk_ele_idx();
166  if (typeid(*isec) == typeid(IntersectionLocal<2,2>)) {
167  //
168  auto il = static_cast<const IntersectionLocal<2,2> *>(isec);
169  ASSERT_EQ( il->size(), 2);
170 
171  arma::vec3 diff = (*il)[0].coords(ele) - (*il)[1].coords(ele);
172  measure_= arma::norm(diff, 2);
173  if (measure_ < 1e-10) return false;
174  } else {
175  // We can save some multiplications moving measure of the master element into
176  // common scale coefficient. However then the meaning of measure_ depends on the case.
177 
178  // multiply by jacobian of transform to real element.
179  double rel_measure = isec->compute_measure();
180  if (rel_measure < 1e-10) return false;
181  measure_= rel_measure * ele.measure() * ele->dim();
182  }
183  return true;
184  }
185 
186 
187  inline uint slave_idx() const
188  {return slave_idx_;}
189 
190  inline double measure() const
191  {return measure_;}
192 
193 private:
195  unsigned int slave_idx_;
196  double measure_;
197 };
198 
199 
200 #endif /* INTERSECTION_HH_ */
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
double measure() const
Computes the measure of the element.
unsigned int dim() const
Definition: elements.h:120
Common base for intersection object.
unsigned int component_ele_idx() const
Returns index of component element.
unsigned int bulk_ele_idx() const
Returns index of bulk element.
virtual double compute_measure() const =0
Class represents intersection of two elements.
bool reinit(const IntersectionLocalBase *isec)
ElementAccessor< 3 > element_accessor(unsigned int idx) const
Create and return ElementAccessor to element of given idx.
Definition: mesh.cc:866
Definition: mesh.h:362
Classes with algorithms for computation of intersections of meshes.
unsigned int uint