← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty-proposed

« back to all changes in this revision

Viewing changes to nfem/lib/geo-interface.cc

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
Import upstream version 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/lib/geo-interface.cc
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
6
/// Rheolef is free software; you can redistribute it and/or modify
 
7
/// it under the terms of the GNU General Public License as published by
 
8
/// the Free Software Foundation; either version 2 of the License, or
 
9
/// (at your option) any later version.
 
10
///
 
11
/// Rheolef is distributed in the hope that it will be useful,
 
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
/// GNU General Public License for more details.
 
15
///
 
16
/// You should have received a copy of the GNU General Public License
 
17
/// along with Rheolef; if not, write to the Free Software
 
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
19
///
 
20
/// =========================================================================
 
21
//
 
22
// geo-interface : provides normal, tangent to subdomains
 
23
//
 
24
// author:
 
25
//      J.Etienne@damtp.cam.ac.uk
 
26
// contribution : splines by B. Balde
 
27
//
 
28
// date:
 
29
//      26/09/06
 
30
//
 
31
#include "georep.h"
 
32
#include "geo-connectivity.h"
 
33
#include "field.h"
 
34
#include <string>
 
35
using namespace std;
 
36
 
 
37
 
 
38
field
 
39
geo::normal (const class space& Nh, const domain& d, const std::string& region) const
 
40
 {
 
41
    interface bc(region);
 
42
    return data().get_normal(Nh, d, false, bc);
 
43
 }
 
44
class field
 
45
geo::normal (const space& Nh, const domain& d, const interface& bc) const
 
46
 {
 
47
    return data().get_normal(Nh, d, false, bc);
 
48
 }
 
49
class field
 
50
geo::tangent (const space& Nh, const domain& d, const string& region) const
 
51
 {
 
52
    interface bc(region);
 
53
    return tangent(Nh, d, bc);
 
54
 }
 
55
class field
 
56
geo::tangent (const space& Nh, const domain& d, const interface& bc) const
 
57
 {
 
58
    field n=data().get_normal(Nh, d, false, bc);
 
59
    field t(n.get_space());
 
60
    t[0]=n[1]; t[1]=-1.*n[0];
 
61
    return t;
 
62
 }
 
63
class field
 
64
geo::axisymmetric_curvature (const space& Nh, const domain& d) const
 
65
 {
 
66
    interface bc;
 
67
    return data().axi_curvature(Nh, d, bc);
 
68
 }
 
69
class field
 
70
geo::axisymmetric_curvature (const space& Nh, const domain& d,
 
71
        const interface& bc) const
 
72
 {
 
73
    return data().axi_curvature(Nh, d, bc);
 
74
 }
 
75
class field
 
76
geo::plane_curvature (const space& Nh, const domain& d,
 
77
        const interface& bc) const
 
78
 {
 
79
    return data().plane_curvature(Nh, d, bc);
 
80
 }
 
81
field
 
82
geo::tangent_spline (const space& Nh, const domain& d, 
 
83
        const interface& bc) const
 
84
 {
 
85
    return data().tangent_spline(Nh, d, bc);
 
86
 }
 
87
class field
 
88
geo::plane_curvature_spline (const space& Nh, const domain& d,
 
89
        const interface& bc) const
 
90
 {
 
91
    return data().plane_curvature_spline(Nh, d, bc);
 
92
 }
 
93
class field
 
94
geo::plane_curvature_quadratic (const space& Nh, const domain& d,
 
95
        const interface& bc) const
 
96
 {
 
97
    return data().plane_curvature_quadratic(Nh, d, bc);
 
98
 }
 
99
 
 
100
void
 
101
georep::set_interface_orientation(const domain& d, const interface& bc) const
 
102
{
 
103
    // set of the vertices of domain d
 
104
    map< size_type, size_type > domain_element;
 
105
    domain::const_iterator iE=d.begin();
 
106
    domain::const_iterator eE=d.end();
 
107
    for (size_type i=0 ; iE!=eE; iE++, i++)
 
108
      { domain_element.insert(pair<size_type,size_type>((*iE).index(),i)); }
 
109
 
 
110
    // find an element of the inside-region with a side on the interface
 
111
    Vector<geo_element>::const_iterator iK, eK;
 
112
    if (bc.region=="")
 
113
     { iK=begin(); eK=end(); }
 
114
    else 
 
115
     { iK=operator[](bc.region).begin(); eK=operator[](bc.region).end(); }
 
116
    map<size_type,size_type>::const_iterator side;
 
117
    map<size_type,size_type>::const_iterator not_found=domain_element.end();
 
118
    bool found=false;
 
119
    for ( ; !found && iK != eK ; iK++ )
 
120
      for (size_type i=0; i<(*iK).n_subgeo((*iK).dimension()-1); i++)
 
121
        {
 
122
          side=domain_element.find((*iK).side(i));
 
123
          if (side!=not_found)
 
124
          {
 
125
            // Get the outward normal to element iK, whose side is on the interface
 
126
            point n = normal ((*iK), i);
 
127
            // remember that this element has this outward normal
 
128
            bc.set_known_normal((*side).second,n);
 
129
            found=true; break;
 
130
          }
 
131
        }
 
132
    check_macro(found, "Interface is not on the boundary of region " << bc.region); 
 
133
}
 
134
 
 
135
// computes a P0 normal (no inverse connectivity reconstructed)
 
136
class field
 
137
georep::get_normal_P0 (const space& Nh, const domain& d, 
 
138
        bool axisymmetric_curvature, const interface& bc) const
 
139
{
 
140
    point er,ez;
 
141
    if (axisymmetric_curvature) 
 
142
     {
 
143
        check_macro(Nh.coordinate_system()=="rz"||Nh.coordinate_system()=="zr",
 
144
                "Non-axisymmetric geometry");
 
145
        if (Nh.coordinate_system()=="rz") { er=point(1,0); ez=point(0,1); }
 
146
        else { er=point(0,1); ez=point(1,0); }
 
147
     }
 
148
    // fill in with zeros in case only part of domain d is bounding region
 
149
    field nh(Nh, Float(0));
 
150
   
 
151
    // set of the vertices of domain d
 
152
    map< size_type, size_type > domain_element;
 
153
    domain::const_iterator iE=d.begin();
 
154
    domain::const_iterator eE=d.end();
 
155
    for (size_type i=0 ; iE!=eE; iE++, i++)
 
156
      { //cerr << (*iE).index() <<":"<<(*iE)[0]<<"-"<<(*iE)[1] << " ";
 
157
        domain_element.insert(pair<size_type,size_type>((*iE).index(),i)); }
 
158
 
 
159
    // go through all elements in region (or whole mesh) and build normal
 
160
    Vector<geo_element>::const_iterator iK, eK;
 
161
    if (bc.region=="")
 
162
     { iK=begin(); eK=end(); }
 
163
    else 
 
164
     { iK=operator[](bc.region).begin(); eK=operator[](bc.region).end(); }
 
165
    tiny_vector<size_type> dofs;
 
166
    map<size_type,size_type>::const_iterator side;
 
167
    map<size_type,size_type>::const_iterator not_found=domain_element.end();
 
168
    bool first_time=true;
 
169
    for ( ; iK != eK ; iK++ )
 
170
      for (size_type i=0; i<(*iK).n_subgeo((*iK).dimension()-1); i++)
 
171
        {
 
172
          side=domain_element.find((*iK).side(i));
 
173
          if (side!=not_found)
 
174
          {
 
175
            // Get the outward normal to element iK, whose side is on the interface
 
176
            // WARNING: straight sided elements assumed, and P0 normal ...!
 
177
            point n = normal ((*iK), i);
 
178
            if (first_time) 
 
179
            { // remember that this element has this outward normal for
 
180
              // use in plane_curvature (gives order on interface)
 
181
              //std::cerr << " $ elenorm: " << _x.begin()[(*iK)[i]] << " || " 
 
182
              //        << _x.begin()[(*iK)[(i+1)%3]] << " || " << _x.begin()[(*iK)[(i+2)%3]]
 
183
              //        << std::endl;
 
184
              bc.set_known_normal((*side).second,n);
 
185
              first_time=false;
 
186
            }
 
187
            if (axisymmetric_curvature) 
 
188
            {
 
189
              // n.er gives the tan of angle of edge to horizontal,
 
190
              // then curvature is calculated as the inverse of the distance
 
191
              // form middle of edge to a point on the axis s.t....
 
192
              point Me=0.5*(_x.begin()[(*iK)[i]]+_x.begin()[(*iK)[(i+1)%3]]); 
 
193
              n= -dot(n,er)/dot(Me,er) *n;
 
194
            }
 
195
            
 
196
            for (size_type id=0; id<Nh.n_component(); id++) 
 
197
            {
 
198
              // get nh dofs corresponding
 
199
              Nh.set_dof((Nh.get_geo().begin())[(*side).second],dofs,id);
 
200
              // set dofs
 
201
              for (size_type idof=0; idof<dofs.size(); idof++)
 
202
                nh.at(dofs[idof])=n[id];
 
203
            }
 
204
          }
 
205
        }
 
206
    return nh;
 
207
}
 
208
// computes a P1 or P2 normal after reordering the elements of the interface (required to be open
 
209
// because of the reordering)
 
210
class field
 
211
georep::get_normal_P1P2 (const space& Nh, const domain& d, const interface& bc) const
 
212
{
 
213
  string approx = Nh.get_approx();
 
214
  // fill in with zeros in case only part of domain d is bounding region
 
215
  field nh(Nh, Float(0));
 
216
  const_iterator_node sommets = begin_node();
 
217
  (*this).sort_interface(d, bc);
 
218
  Vector<size_t>::const_iterator bI = bc.begin();
 
219
  Vector<size_t>::const_iterator eI = bc.end();
 
220
  //we take care of the first two edges
 
221
  geo_element e0 = d.at(*bI++);
 
222
  geo_element e1 = d.at(*bI++);
 
223
  //the edges are sorted, but not necessarily the nodes
 
224
  size_type OrderedLocalIndex0 [2];
 
225
  size_type OrderedLocalIndex1 [2];
 
226
  if (e0[0]==e1[1]) {
 
227
    OrderedLocalIndex0[0] = 1;
 
228
    OrderedLocalIndex0[1] = 0;
 
229
    OrderedLocalIndex1[0] = 1;
 
230
    OrderedLocalIndex1[1] = 0;
 
231
  } 
 
232
  else if (e0[1]==e1[1]) {
 
233
    OrderedLocalIndex0[0] = 0;
 
234
    OrderedLocalIndex0[1] = 1;
 
235
    OrderedLocalIndex1[0] = 1;
 
236
    OrderedLocalIndex1[1] = 0;
 
237
  }
 
238
  else if (e0[1]==e1[0]) {
 
239
    OrderedLocalIndex0[0] = 0;
 
240
    OrderedLocalIndex0[1] = 1;
 
241
    OrderedLocalIndex1[0] = 0;
 
242
    OrderedLocalIndex1[1] = 1;
 
243
  }
 
244
  else /* e0[0]==e1[0] */ {
 
245
    OrderedLocalIndex0[0] = 1;
 
246
    OrderedLocalIndex0[1] = 0;
 
247
    OrderedLocalIndex1[0] = 0;
 
248
    OrderedLocalIndex1[1] = 1;
 
249
  }
 
250
  Float coord0[2][2]; 
 
251
  Float coord1[2][2]; 
 
252
  for (size_type i = 0 ; i < 2 ; i++) {
 
253
    for (size_type j = 0 ; j < 2 ; j++) {
 
254
      coord0[i][j] = sommets[e0[OrderedLocalIndex0[i]]][j];
 
255
      coord1[i][j] = sommets[e1[OrderedLocalIndex1[i]]][j];
 
256
    }
 
257
  }
 
258
  point n0(-(coord0[1][1]-coord0[0][1]), coord0[1][0]-coord0[0][0]);
 
259
  point n1(-(coord1[1][1]-coord1[0][1]), coord1[1][0]-coord1[0][0]);
 
260
  point n_avrg = n0 + n1;
 
261
  n_avrg = n_avrg/norm(n_avrg);
 
262
  Float edge0 = norm(n0);
 
263
  Float edge1 = norm(n1);
 
264
  tiny_vector<size_t> dof;
 
265
  size_type d0,d1;
 
266
  for (size_type id = 0 ; id < Nh.n_component(); id++){
 
267
    Nh.set_global_dof(e0, dof, id);
 
268
    d0 = dof[OrderedLocalIndex0[0]];
 
269
    d1 = dof[OrderedLocalIndex0[1]];
 
270
    if (bc.first_normal_imposed()||bc.first_tangent_imposed())
 
271
      nh.at(Nh.domain_index(d0)) = (bc.first_normal())(id);
 
272
    else {
 
273
      nh.at(Nh.domain_index(d0)) = n0(id)/edge0;
 
274
      //warning_macro("first normal is not imposed. Use set_first_normal() if you want to impose it");
 
275
    }
 
276
    if (approx == "P2")
 
277
      nh.at(Nh.domain_index(dof[2])) = n0(id)/edge0;
 
278
    nh.at(Nh.domain_index(d1)) = n_avrg(id);
 
279
    if (approx == "P2") {
 
280
      Nh.set_global_dof(e1, dof, id);
 
281
      nh.at(Nh.domain_index(dof[2])) = n1(id)/edge1;
 
282
    }
 
283
  }
 
284
  n0 = n1;
 
285
  size_type indexp = e1[OrderedLocalIndex1[1]];
 
286
  for( ; bI != eI ; bI++) {
 
287
    geo_element e = d.at(*bI);
 
288
    if (indexp == e[0]) {
 
289
      OrderedLocalIndex1[0] = 0;
 
290
      OrderedLocalIndex1[1] = 1;
 
291
    } else {
 
292
      OrderedLocalIndex1[0] = 1;
 
293
      OrderedLocalIndex1[1] = 0;
 
294
    }
 
295
    for (size_type i = 0 ; i < 2 ; i++) {
 
296
      for (size_type j = 0 ; j < 2 ; j++) {
 
297
        coord1[i][j] = sommets[e[OrderedLocalIndex1[i]]][j];
 
298
      }
 
299
    }
 
300
    n1 = point(-(coord1[1][1]-coord1[0][1]), (coord1[1][0]-coord1[0][0]));
 
301
    n_avrg = n0 + n1;
 
302
    edge0 = norm(n0);
 
303
    edge1 = norm(n1);
 
304
    n_avrg = n_avrg/sqrt(dot(n_avrg,n_avrg));
 
305
    for (size_type id = 0 ; id < Nh.n_component(); id++){
 
306
      Nh.set_global_dof(e, dof, id);
 
307
      d0 = dof[OrderedLocalIndex1[0]];
 
308
      d1 = dof[OrderedLocalIndex1[1]];
 
309
      nh.at(Nh.domain_index(d0)) = n_avrg(id);
 
310
      if (approx == "P2")
 
311
        nh.at(Nh.domain_index(dof[2])) = n1(id)/edge1;
 
312
    }
 
313
    indexp = e[OrderedLocalIndex1[1]];
 
314
    n0 = n1;
 
315
  }
 
316
  geo_element e = d.at(*(--bI));
 
317
  for (size_type id = 0 ; id < Nh.n_component(); id++){
 
318
    Nh.set_global_dof(e, dof, id);
 
319
    d1 = dof[OrderedLocalIndex1[1]];
 
320
    if (bc.last_normal_imposed()||bc.last_tangent_imposed())
 
321
      nh.at(Nh.domain_index(d1)) = (bc.last_normal())(id);
 
322
    else {
 
323
      nh.at(Nh.domain_index(d1)) = n1(id)/edge1;
 
324
      //warning_macro("last normal is not imposed. Use set_last_normal() if you want to impose it");
 
325
    }
 
326
  }
 
327
  return nh;
 
328
}
 
329
// computes a P0, P1, or P2 normal 
 
330
class field
 
331
georep::get_normal (const space& Nh, const domain& d, 
 
332
        bool axisymmetric_curvature, const interface& bc) const
 
333
{
 
334
  string approx = Nh.get_approx();
 
335
  check_macro((approx=="P0"||approx=="P1"||approx=="P2"), 
 
336
              approx + " normal not implemented.");
 
337
  if (approx=="P0") {
 
338
    return get_normal_P0 (Nh, d, axisymmetric_curvature, bc);
 
339
  }
 
340
  check_macro(!axisymmetric_curvature, "axisymetric curvature not implemented for " + approx + " normals.");
 
341
  return get_normal_P1P2 (Nh, d, bc);
 
342
}
 
343
inline
 
344
point ortho2d(const point& AB)
 
345
{
 
346
    return point(-AB[1],AB[0]);
 
347
}
 
348
inline
 
349
point normal2d(const point& AB)
 
350
{
 
351
    return point(-AB[1],AB[0])/norm(AB);
 
352
}
 
353
 
 
354
inline
 
355
Float curvature2d(const point& AB, const point& BC)
 
356
// calculates the curvature at B using a quadratic approximation
 
357
// with sA=-|AB|, sB=0, sC=|BC|, we get
 
358
//  x'y''-x''y' = -2 |AB x BC| / ( sA sC (sC-sA) )
 
359
// (which is indep of s), and
 
360
//  (x'^2+y'^2)(s=0) = |AC|^2 /(sC-sA)^2
 
361
 {
 
362
   Float sA=-norm(AB); 
 
363
   Float sC=norm(BC);
 
364
   return -2*(sC-sA)*vect2d(AB,BC) / sA / sC / norm2(AB+BC);
 
365
 }
 
366
inline
 
367
point vcurvature2d(const point& AB, const point& BC)
 
368
// calculates the curvature at B using a quadratic approximation
 
369
 {
 
370
   Float sA2=norm2(AB); Float sA=-sqrt(sA2); 
 
371
   Float sC2=norm2(BC); Float sC=sqrt(sC2); 
 
372
   return 2*(sC-sA)*vect2d(AB,BC) / sA2 / sC2 / pow(norm2(AB+BC), Float(1.5))
 
373
                *(sA2*ortho2d(BC)+sC2*ortho2d(AB));
 
374
 }
 
375
inline
 
376
point normal2d(const point& AB, const point& BC)
 
377
// calculates the normal at B using a quadratic approximation
 
378
 {
 
379
   Float sA2=norm2(AB); Float sA=-sqrt(sA2); 
 
380
   Float sC2=norm2(BC); Float sC=sqrt(sC2); 
 
381
   return -1 / sA / sC / norm(AB+BC) *(sA2*ortho2d(BC)+sC2*ortho2d(AB));
 
382
 }
 
383
inline
 
384
Float curvature2d_axi(const point& M, const point& AB, const point& er)
 
385
// n.er = |AB x er|/|AB| gives the tan of angle of edge to horizontal,
 
386
// then curvature is calculated as the inverse of the distance
 
387
// form middle of edge to a point on the axis s.t....
 
388
 {
 
389
   return -vect2d(AB,er)/dot(M,er)/norm(AB);
 
390
 }
 
391
 
 
392
/*
 
393
// calculates the inverse of the radius of a circle going through A, B and C
 
394
// curvature sign is determined independently
 
395
{
 
396
    if (abs(vect2d(AB,BC))<sqrt(numeric_limits<Float>::epsilon()))
 
397
        // points are aligned
 
398
        return 0;
 
399
    Float sign=1;
 
400
    if (vect2d(AB,BC)<0) sign=-1;
 
401
    // let M be the middle of BC:
 
402
    Float dist_OM=-.5*dot(AB+BC,AB)/dot(normal2d(BC),AB);
 
403
    return sign/sqrt(sqr(dist_OM)+norm2(BC)/4.);
 
404
}
 
405
inline
 
406
point vcurvature2d(const point& AB, const point& BC)
 
407
{
 
408
    if (abs(vect2d(AB,BC))<sqrt(numeric_limits<Float>::epsilon()))
 
409
        // points are aligned
 
410
        return point(0);
 
411
    Float signed_dist_OM=-.5*dot(AB+BC,AB)/dot(normal2d(BC),AB);
 
412
    return .5*BC+signed_dist_OM*normal2d(BC);
 
413
}
 
414
*/
 
415
 
 
416
void
 
417
georep::sort_interface (const domain& d, const interface& bc) const
 
418
 {
 
419
    if (bc.is_sorted()) return;
 
420
    // Initialize interface boundary condition 
 
421
    if (!bc.initialized())
 
422
        set_interface_orientation(d, bc);
 
423
    // Sort contiguous elements of the interface
 
424
    //  first, group them in contiguous stubs listing element #s that are merged
 
425
    //  together as the jigsaw progresses...
 
426
    typedef list< size_type > stub;
 
427
    typedef list< stub > set_of_stubs;
 
428
    set_of_stubs S;
 
429
    set_of_stubs::iterator iS;
 
430
    set_of_stubs::iterator eS;
 
431
    set_of_stubs::iterator iS_attached;
 
432
    // lists of the first and last element in stubs -- misleadingly of type 'stub' as
 
433
    // well, these contain *vertex indices*, not *element* ones. 
 
434
    stub F, L; 
 
435
    stub::iterator iF, iL, iF_attached, iL_attached;
 
436
    domain::const_iterator bd=d.begin();
 
437
    domain::const_iterator ed=d.end();
 
438
    size_type d_index=0;
 
439
    bool actual_order_is_set=false;
 
440
    for (domain::const_iterator id=bd; id!=ed; id++, d_index++)
 
441
    {
 
442
        bool attached=false;
 
443
        bool over=false;
 
444
        int attached_dir=-1;
 
445
        for (size_type i=0; i<2; i++) // examine both orders for attachment
 
446
        for (iS=S.begin(), eS=S.end(),
 
447
             iF=F.begin(),
 
448
             iL=L.begin();
 
449
             !over && iS!=eS; iS++, iF++, iL++)
 
450
        {
 
451
            if (*iF==(*id)[i]) 
 
452
            {
 
453
            // first element of the stub has a common vertex with *id, then
 
454
            // attach *id to front of stub *iS
 
455
                if (!attached) 
 
456
                // other end of id is still free
 
457
                {
 
458
                  (*iS).push_front(d_index);
 
459
                  *iF=((*id)[(i+1)%2]);
 
460
                  attached=true;
 
461
                  iS_attached=iS;
 
462
                  iF_attached=iF;
 
463
                  iL_attached=iL;
 
464
                  attached_dir=i;
 
465
                }
 
466
                else
 
467
                if (iS_attached!=iS) 
 
468
                // other end is already attached to a stub: merge them
 
469
                {
 
470
                  if (i==0) // the other stub is in the other order
 
471
                  {
 
472
                    *iF=*iL_attached;
 
473
                    (*iS_attached).reverse();
 
474
                  }
 
475
                  else *iF=*iF_attached;   
 
476
                  (*iS).splice((*iS).begin(),*iS_attached);
 
477
                  (*iS_attached).~stub();
 
478
                  F.erase(iF_attached);
 
479
                  L.erase(iL_attached);
 
480
                  S.erase(iS_attached);
 
481
                  over=true;
 
482
                }
 
483
            }
 
484
//          if (bd[*(*iS).rbegin()][1]==(*id)[i])
 
485
            if (*iL==(*id)[i])
 
486
            {
 
487
            // last element of the stub has a common vertex with *id, then
 
488
            // attach *id to back of stub *iS 
 
489
                if (!attached)
 
490
                {
 
491
                  (*iS).push_back(d_index);
 
492
                  *iL=((*id)[(i+1)%2]);
 
493
                  attached=true;
 
494
                  iS_attached=iS;
 
495
                  iF_attached=iF;
 
496
                  iL_attached=iL;
 
497
                  attached_dir=i;
 
498
                }
 
499
                else
 
500
                if (iS_attached!=iS) 
 
501
                {
 
502
                  if (i==1) // the other stub is in the other order
 
503
                  {
 
504
                    *iL=*iF_attached;
 
505
                    (*iS_attached).reverse();
 
506
                  }
 
507
                  else *iL=*iL_attached;
 
508
                  (*iS_attached).splice((*iS_attached).begin(), *iS);
 
509
                  (*iS).~stub();
 
510
                  F.erase(iF_attached);
 
511
                  L.erase(iL_attached);
 
512
                  S.erase(iS);
 
513
                  over=true;
 
514
                }
 
515
            }
 
516
        }
 
517
        if (!attached) // we have to create a new stub for this element
 
518
        {
 
519
            stub* newstub =new stub;
 
520
            (*newstub).push_front(d_index);
 
521
            S.push_front(*newstub);
 
522
            F.push_front((*id)[0]);
 
523
            L.push_front((*id)[1]);
 
524
            attached_dir=0;
 
525
        }
 
526
        if (d_index==bc.element_of_known_normal)
 
527
        {
 
528
            check_macro(attached_dir!=-1,"Attached direction not set !" << S.size());
 
529
            bc.set_actual_order( ((attached_dir==0)?1:-1)
 
530
                *(_x.begin()[(*id)[1]] - _x.begin()[(*id)[0]]));
 
531
            //      cerr << "$$ Actual order : " << ((attached_dir==0)?1:-1)*(_x.begin()[(*id)[1]] - _x.begin()[(*id)[0]]) << "( dir " << ((attached_dir==0)?1:-1) << " )"<< endl;
 
532
            actual_order_is_set=true;
 
533
        }
 
534
    }
 
535
    check_macro (actual_order_is_set, "Element whose normal is marked " 
 
536
        << bc.element_of_known_normal << " not found!");
 
537
    check_macro (S.size()==1, "The interface has to be a connected set!");
 
538
    check_macro ((*F.begin())!=(*L.begin()), "Closed interfaces not supported yet!");
 
539
    check_macro ((*S.begin()).size()==d.size(), "Not all elements found?");
 
540
    if (bc.is_reversed())
 
541
      bc.set_sorted( (*S.begin()).size(), (*S.begin()).rbegin(), (*S.begin()).rend());
 
542
    else
 
543
      bc.set_sorted( (*S.begin()).size(), (*S.begin()).begin(), (*S.begin()).end());
 
544
    
 
545
    // get global dof number for first and last dof
 
546
    Vector<size_t>::const_iterator iAB=bc.begin();
 
547
    Vector<size_t>::const_iterator iBC=bc.begin(); iBC++;
 
548
    size_t AB0=bd[*iAB][0];
 
549
    size_t AB1=bd[*iAB][1];
 
550
    size_t BC0=bd[*iBC][0];
 
551
    size_t BC1=bd[*iBC][1];
 
552
    // check that nodes are in the same order as elements, otherwise reverse
 
553
    if (AB1==BC0||AB1==BC1) { bc.set_first_node(*iAB, 0); }
 
554
    else
 
555
    if (AB0==BC0||AB0==BC1) { bc.set_first_node(*iAB, 1); }
 
556
    else error_macro("Fatal: non consecutive edges!");
 
557
 
 
558
    iAB=bc.end(); iAB--;
 
559
    iBC=iAB; iBC--;
 
560
    AB0=bd[*iAB][0];
 
561
    AB1=bd[*iAB][1];
 
562
    BC0=bd[*iBC][0];
 
563
    BC1=bd[*iBC][1];
 
564
    // check that nodes are in the same order as elements, otherwise reverse
 
565
    if (AB1==BC0||AB1==BC1) { bc.set_last_node(*iAB,0); }
 
566
    else
 
567
    if (AB0==BC0||AB0==BC1) { bc.set_last_node(*iAB,1); }
 
568
    else error_macro("Fatal: non consecutive edges!");
 
569
 }
 
570
field
 
571
georep::plane_curvature (const space& Nh, const domain& d, 
 
572
        const interface& bc) const
 
573
{
 
574
    return plane_curvature_spline (Nh, d, bc);
 
575
}
 
576
 
 
577
field
 
578
georep::tangent_spline (const space& Nh, const domain& d, 
 
579
        const interface& bc) const
 
580
{
 
581
  check_macro(Nh.get_geo() .dimension()==2,"plane curvature is for 2D geometries only");
 
582
  check_macro(Nh.n_component ()==2,"2-component vector space needed");
 
583
  // Sort contiguous elements of the interface
 
584
  // first ,group them in contiguous stubs listing element #s that are merged
 
585
  // together as the jidsaw progresses....
 
586
  sort_interface(d,bc);
 
587
  // get coodinate of lagrange interpolation points on edge elements :
 
588
  vector<point> hat_node;
 
589
  domain::const_iterator bd=d.begin () ;
 
590
  Nh.get_basis().hat_node(*bd,hat_node) ;
 
591
 
 
592
  Vector<size_t>::const_iterator iAB=bc.begin() ;
 
593
  Vector<size_t>::const_iterator iBC=bc.begin() ;  iBC++;
 
594
  Vector<size_t>::const_iterator the_end=bc.end() ;
 
595
 
 
596
  // calculate interface curvature using a parametric spline
 
597
  
 
598
  field kh(Nh,float(0));
 
599
  tiny_vector<size_type> dofs;
 
600
 
 
601
  point AB0=_x.begin()[bd[*iAB][0]];
 
602
  point AB1=_x.begin()[bd[*iAB][1]];
 
603
  point BC0=_x.begin()[bd[*iBC][0]];
 
604
  point BC1=_x.begin()[bd[*iBC][1]];
 
605
  point AB,BC;
 
606
  bool rev_AB;
 
607
 
 
608
  /* calulate zi
 
609
     compute the hi and the bi */
 
610
  size_type n=d.size();
 
611
  vector<Float> u(n+1),h(n);
 
612
  vector<point> v(n+1),b(n+1),z(n+1) ,y(n+1),w(n+1);
 
613
  size_type i;
 
614
  for (i=0 ; iBC!=the_end; iAB++,iBC++,i++)
 
615
    {
 
616
        // check that nodes are in the same order as elements ,otherwise reverse 
 
617
        BC0=_x.begin()[bd[*iBC][0]];
 
618
        BC1=_x.begin()[bd[*iBC][1]];
 
619
        if(AB1==BC0) { rev_AB=false; AB=AB1-AB0 ; BC=BC1-BC0; }
 
620
        else
 
621
        if(AB0==BC0) { rev_AB=true; AB=AB0-AB1; BC=BC1-BC0 ; }
 
622
        else
 
623
        if(AB1==BC1) { rev_AB=false;AB=AB1-AB0; BC=BC0-BC1 ; }
 
624
        else
 
625
        {
 
626
                check_macro (AB0==BC1, "Fatal :non consecutive edges!");
 
627
                rev_AB=true; AB=AB0-AB1; BC=BC0-BC1;
 
628
        }
 
629
        h[i]=norm(AB);
 
630
        b[i]=1./h[i]*AB;
 
631
        AB0=BC0; AB1=BC1 ;AB=BC;
 
632
         
 
633
    }
 
634
   
 
635
   h[n-1]=norm(BC); b[n-1]=1./h[n-1]*BC;
 
636
  // GAUSSIAN ELIMINATION 
 
637
  size_type n_loop=n;
 
638
  size_type start_loop=1;
 
639
  point t0;
 
640
  point tn;
 
641
  if (bc.first_tangent_imposed())
 
642
    {
 
643
     t0=bc.first_tangent();
 
644
     u[0]=2*(h[0]);
 
645
     v[0]=6*(b[0]-t0);
 
646
     u[1]=2*(h[0]+h[1])-sqr(h[0])/u[0];
 
647
     v[1]=6*(b[1]-b[0])- (h[0]/u[0])*v[0];
 
648
     start_loop =0;
 
649
    }
 
650
  else 
 
651
    {
 
652
     u[1]=2*(h[0]+h[1]);
 
653
     v[1]=6*(b[1]-b[0]);
 
654
     start_loop=1 ;
 
655
    }
 
656
  for (i=2; i<=(n-1); i++)
 
657
     { 
 
658
       u[i]=2.*(h[i-1]+h[i])-sqr(h[i-1])/u[i-1];
 
659
       v[i]=6.*(b[i]-b[i-1])-h[i-1]/u[i-1]*v[i-1];
 
660
     }
 
661
  if (bc.last_tangent_imposed())
 
662
   { 
 
663
    tn=bc.last_tangent();
 
664
    n_loop=n ;
 
665
    u[n]=2*h[n-1]-sqr(h[n-1])/u[n-1] ;
 
666
    v[n]= 6*(tn-b[n-1])-h[n-1]/u[n-1]*v[n-1];
 
667
   }
 
668
  else 
 
669
   { 
 
670
    n_loop=n-1;
 
671
    u[n]=2*h[n-1]-sqr(h[n-1])/u[n-1] ;
 
672
    v[n]= 6*(tn-b[n-1])-h[n-1]/u[n-1]*v[n-1];
 
673
   }
 
674
      
 
675
 
 
676
  //back substitution 
 
677
  if (!bc.last_tangent_imposed()) {z[n]=point(0,0);}
 
678
  else z[n]=1./u[n]*v[n];
 
679
     for (int j=n-1; j>=int(start_loop); j--)
 
680
       {
 
681
       z[j]=1./u[j]*( v[j] -h[j]*z[j+1] );
 
682
     
 
683
       w[j]=-h[j]/6 *(z[j+1]-z[j]) +b[j];
 
684
      
 
685
       }
 
686
  if (!bc.first_tangent_imposed()){z[0]=point(0,0);w[0]=-h[0]/6 *(z[1]-z[0]) +b[0];}
 
687
   
 
688
  // curvature calculation
 
689
  iAB=bc.begin() ;
 
690
  iBC=bc.begin() ;  iBC++;
 
691
  AB0=_x.begin()[bd[*iAB][0]];
 
692
  AB1=_x.begin()[bd[*iAB][1]];
 
693
 
 
694
  point tang;
 
695
  for (i=0 ; iBC!=the_end; iAB++,iBC++,i++)
 
696
   {
 
697
        BC0=_x.begin()[bd[*iBC][0]];
 
698
        BC1=_x.begin()[bd[*iBC][1]];
 
699
        // check that nodes are in the same order as elements ,otherwise reverse 
 
700
 
 
701
        if(AB1==BC0) { rev_AB=false; AB=AB1-AB0 ; BC=BC1-BC0; }
 
702
        else
 
703
        if(AB0==BC0) { rev_AB=true; AB=AB0-AB1; BC=BC1-BC0 ; }
 
704
        else
 
705
        if(AB1==BC1) { rev_AB=false;AB=AB1-AB0; BC=BC0-BC1 ; }
 
706
        else
 
707
        {
 
708
                check_macro (AB0==BC1, "Fatal :non consecutive edges!");
 
709
                rev_AB=true; AB=AB0-AB1; BC=BC0-BC1;
 
710
        }
 
711
        for (size_type j=0;j<Nh.n_component();j++)
 
712
          {
 
713
                // get nh dofs corresponding
 
714
 
 
715
                Nh.set_dof((Nh.get_geo().begin())[*iAB],dofs,j);
 
716
                //set dofs
 
717
 
 
718
                for (size_type idof=0;idof<dofs.size();idof++)
 
719
                  {
 
720
                   if(!rev_AB)
 
721
                     {
 
722
                        tang=-sqr(1-hat_node[idof][0])*z[i]/2*h[i]  
 
723
                                    +sqr(hat_node[idof][0])*z[i+1]/2*h[i] 
 
724
                                    +w[i];
 
725
                     }
 
726
                   else
 
727
                     {
 
728
                        //cerr << (1-hat_node[idof][0]) << "*" << z[i][j] 
 
729
                            // << " + " << hat_node[idof][0] << "*" << z[i+1][j] << endl;
 
730
                        tang=sqr(1-hat_node[idof][0])* z[i+1]/2*h[i] 
 
731
                                   -sqr(hat_node[idof][0])*z[i]/2*h[i] 
 
732
                                   +w[i];
 
733
                     }
 
734
                   kh.at(dofs[idof])= tang[j];
 
735
                 
 
736
                 } // idof
 
737
          } // j<n_compo
 
738
        AB0=BC0; AB1=BC1 ;AB=BC;
 
739
  } // iAB
 
740
 // don't need to check that nodes are in the same order as elements,as BC was
 
741
 
 
742
 rev_AB=(AB==AB0 -AB1);
 
743
 i=n-1;
 
744
 for (size_type j=0;j<Nh.n_component();j++)
 
745
 {
 
746
                // get nh dofs corresponding
 
747
 
 
748
                Nh.set_dof((Nh.get_geo().begin())[*iAB],dofs,j);
 
749
                //set dofs
 
750
              
 
751
                for (size_type idof=0;idof<dofs.size();idof++)
 
752
                 {
 
753
                
 
754
                   if(!rev_AB)
 
755
                     {
 
756
                        tang=-sqr(1-hat_node[idof][0])*z[i]/2*h[i] 
 
757
                                    +sqr(hat_node[idof][0])* z[i+1]/2*h[i] 
 
758
                                    +w[i];
 
759
                     }
 
760
                   else
 
761
                     {
 
762
                        //cerr << (1-hat_node[idof][0]) << "*" << z[i][j] 
 
763
                            // << " + " << hat_node[idof][0] << "*" << z[i+1][j] << endl;
 
764
                        tang=sqr(1-hat_node[idof][0])*z[i+1]/2*h[i] 
 
765
                                   -sqr(hat_node[idof][0])*z[i]/2*h[i] 
 
766
                                   +w[i];
 
767
                     }
 
768
                   kh.at(dofs[idof])= tang[j];
 
769
                  
 
770
                 }
 
771
  }
 
772
  return kh;
 
773
}
 
774
 
 
775
field
 
776
georep::plane_curvature_spline (const space& Nh, const domain& d, 
 
777
        const interface& bc) const
 
778
{
 
779
  check_macro(Nh.get_geo() .dimension()==2,"plane curvature is for 2D geometries only");
 
780
  check_macro(Nh.n_component ()==2,"2-component vector space needed");
 
781
  // Sort contiguous elements of the interface
 
782
  // first ,group them in contiguous stubs listing element #s that are merged
 
783
  // together as the jidsaw progresses....
 
784
  sort_interface(d,bc);
 
785
  // get coodinate of lagrange interpolation points on edge elements :
 
786
  vector<point> hat_node;
 
787
  domain::const_iterator bd=d.begin () ;
 
788
  Nh.get_basis().hat_node(*bd,hat_node) ;
 
789
 
 
790
  Vector<size_t>::const_iterator iAB=bc.begin() ;
 
791
  Vector<size_t>::const_iterator iBC=bc.begin() ;  iBC++;
 
792
  Vector<size_t>::const_iterator the_end=bc.end() ;
 
793
 
 
794
  // calculate interface curvature using a parametric spline
 
795
  
 
796
  field kh(Nh,float(0));
 
797
  tiny_vector<size_type> dofs;
 
798
 
 
799
  point AB0=_x.begin()[bd[*iAB][0]];
 
800
  point AB1=_x.begin()[bd[*iAB][1]];
 
801
  point BC0=_x.begin()[bd[*iBC][0]];
 
802
  point BC1=_x.begin()[bd[*iBC][1]];
 
803
  point AB,BC;
 
804
  bool rev_AB;
 
805
 
 
806
  /* calulate zi
 
807
     compute the hi and the bi */
 
808
  size_type n=d.size();
 
809
  vector<Float> u(n+1),h(n);
 
810
  vector<point> v(n+1),b(n+1),z(n+1) ,y(n+1),w(n+1);
 
811
  size_type i;
 
812
  for (i=0 ; iBC!=the_end; iAB++,iBC++,i++)
 
813
    {
 
814
        // check that nodes are in the same order as elements ,otherwise reverse 
 
815
        BC0=_x.begin()[bd[*iBC][0]];
 
816
        BC1=_x.begin()[bd[*iBC][1]];
 
817
        if(AB1==BC0) { rev_AB=false; AB=AB1-AB0 ; BC=BC1-BC0; }
 
818
        else
 
819
        if(AB0==BC0) { rev_AB=true; AB=AB0-AB1; BC=BC1-BC0 ; }
 
820
        else
 
821
        if(AB1==BC1) { rev_AB=false;AB=AB1-AB0; BC=BC0-BC1 ; }
 
822
        else
 
823
        {
 
824
                check_macro (AB0==BC1, "Fatal :non consecutive edges!");
 
825
                rev_AB=true; AB=AB0-AB1; BC=BC0-BC1;
 
826
        }
 
827
        h[i]=norm(AB);
 
828
        b[i]=1./h[i]*AB;
 
829
        AB0=BC0; AB1=BC1 ;AB=BC;
 
830
         
 
831
    }
 
832
   
 
833
   h[n-1]=norm(BC); b[n-1]=1./h[n-1]*BC;
 
834
  // GAUSSIAN ELIMINATION 
 
835
  size_type n_loop=n;
 
836
  size_type start_loop=1;
 
837
  point t0;
 
838
  point tn;
 
839
  if (bc.first_tangent_imposed())
 
840
    {
 
841
     t0=bc.first_tangent();
 
842
     u[0]=2*(h[0]);
 
843
     v[0]=6*(b[0]-t0);
 
844
     u[1]=2*(h[0]+h[1])-sqr(h[0])/u[0];
 
845
     v[1]=6*(b[1]-b[0])- (h[0]/u[0])*v[0];
 
846
     start_loop =0;
 
847
    }
 
848
  else 
 
849
    {
 
850
     u[1]=2*(h[0]+h[1]);
 
851
     v[1]=6*(b[1]-b[0]);
 
852
     start_loop=1 ;
 
853
    }
 
854
  for (i=2; i<=(n-1); i++)
 
855
     { 
 
856
       u[i]=2.*(h[i-1]+h[i])-sqr(h[i-1])/u[i-1];
 
857
       v[i]=6.*(b[i]-b[i-1])-h[i-1]/u[i-1]*v[i-1];
 
858
     }
 
859
  if (bc.last_tangent_imposed())
 
860
   { 
 
861
    tn=bc.last_tangent();
 
862
    n_loop=n ;
 
863
    u[n]=2*h[n-1]-sqr(h[n-1])/u[n-1] ;
 
864
    v[n]= 6*(tn-b[n-1])-h[n-1]/u[n-1]*v[n-1];
 
865
   }
 
866
  else 
 
867
   { 
 
868
    n_loop=n-1;
 
869
    u[n]=2*h[n-1]-sqr(h[n-1])/u[n-1] ;
 
870
    v[n]= 6*(tn-b[n-1])-h[n-1]/u[n-1]*v[n-1];
 
871
   }
 
872
      
 
873
 
 
874
  //back substitution 
 
875
  if (!bc.last_tangent_imposed()) {z[n]=point(0,0);}
 
876
  else z[n]=1./u[n]*v[n];
 
877
     for (int j=n-1; j>=int(start_loop); j--)
 
878
       {
 
879
       z[j]=1./u[j]*( v[j] -h[j]*z[j+1] );
 
880
     
 
881
       w[j]=-h[j]/6 *(z[j+1]-z[j]) +b[j];
 
882
      
 
883
       }
 
884
  if (!bc.first_tangent_imposed()){z[0]=point(0,0);w[0]=-h[0]/6 *(z[1]-z[0]) +b[0];}
 
885
   
 
886
  // curvature calculation
 
887
  iAB=bc.begin() ;
 
888
  iBC=bc.begin() ;  iBC++;
 
889
  AB0=_x.begin()[bd[*iAB][0]];
 
890
  AB1=_x.begin()[bd[*iAB][1]];
 
891
 
 
892
  point z_at_dof, devy_at_dof, n1;
 
893
  for (i=0 ; iBC!=the_end; iAB++,iBC++,i++)
 
894
   {
 
895
        BC0=_x.begin()[bd[*iBC][0]];
 
896
        BC1=_x.begin()[bd[*iBC][1]];
 
897
        // check that nodes are in the same order as elements ,otherwise reverse 
 
898
 
 
899
        if(AB1==BC0) { rev_AB=false; AB=AB1-AB0 ; BC=BC1-BC0; }
 
900
        else
 
901
        if(AB0==BC0) { rev_AB=true; AB=AB0-AB1; BC=BC1-BC0 ; }
 
902
        else
 
903
        if(AB1==BC1) { rev_AB=false;AB=AB1-AB0; BC=BC0-BC1 ; }
 
904
        else
 
905
        {
 
906
                check_macro (AB0==BC1, "Fatal :non consecutive edges!");
 
907
                rev_AB=true; AB=AB0-AB1; BC=BC0-BC1;
 
908
        }
 
909
        for (size_type j=0;j<Nh.n_component();j++)
 
910
          {
 
911
                // get nh dofs corresponding
 
912
 
 
913
                Nh.set_dof((Nh.get_geo().begin())[*iAB],dofs,j);
 
914
                //set dofs
 
915
 
 
916
                for (size_type idof=0;idof<dofs.size();idof++)
 
917
                  {
 
918
                   if(!rev_AB)
 
919
                     {
 
920
                        z_at_dof= (1-hat_node[idof][0])* z[i] + hat_node[idof][0]*z[i+1] ;
 
921
                        devy_at_dof=-sqr(1-hat_node[idof][0])*z[i]/2*h[i]  
 
922
                                    +sqr(hat_node[idof][0])*z[i+1]/2*h[i] 
 
923
                                    +w[i];
 
924
                        n1=point(devy_at_dof[1],-devy_at_dof[0]);
 
925
                     }
 
926
                   else
 
927
                     {
 
928
                        //cerr << (1-hat_node[idof][0]) << "*" << z[i][j] 
 
929
                            // << " + " << hat_node[idof][0] << "*" << z[i+1][j] << endl;
 
930
                        z_at_dof=(1-hat_node[idof][0])* z[i+1] + hat_node[idof][0] *z[i];
 
931
                        devy_at_dof=sqr(1-hat_node[idof][0])* z[i+1]/2*h[i] 
 
932
                                   -sqr(hat_node[idof][0])*z[i]/2*h[i] 
 
933
                                   +w[i];
 
934
                        n1=point(devy_at_dof[1],-devy_at_dof[0]);
 
935
                     }
 
936
                   kh.at(dofs[idof])= -1./(sqr(sqr(devy_at_dof[0])+sqr(devy_at_dof[1])))
 
937
                                *(z_at_dof[1]*devy_at_dof[0]-z_at_dof[0]*devy_at_dof[1])
 
938
                                *n1[j];
 
939
                 
 
940
                 } // idof
 
941
          } // j<n_compo
 
942
        AB0=BC0; AB1=BC1 ;AB=BC;
 
943
  } // iAB
 
944
 // don't need to check that nodes are in the same order as elements,as BC was
 
945
 
 
946
 rev_AB=(AB==AB0 -AB1);
 
947
 i=n-1;
 
948
 for (size_type j=0;j<Nh.n_component();j++)
 
949
 {
 
950
                // get nh dofs corresponding
 
951
 
 
952
                Nh.set_dof((Nh.get_geo().begin())[*iAB],dofs,j);
 
953
                //set dofs
 
954
              
 
955
                for (size_type idof=0;idof<dofs.size();idof++)
 
956
                 {
 
957
                
 
958
                   if(!rev_AB)
 
959
                     {
 
960
                        z_at_dof= (1-hat_node[idof][0])* z[i] + hat_node[idof][0] *z[i+1];
 
961
                        devy_at_dof=-sqr(1-hat_node[idof][0])*z[i]/2*h[i] 
 
962
                                    +sqr(hat_node[idof][0])* z[i+1]/2*h[i] 
 
963
                                    +w[i];
 
964
                        n1=point(devy_at_dof[1],-devy_at_dof[0]);
 
965
                     }
 
966
                   else
 
967
                     {
 
968
                        //cerr << (1-hat_node[idof][0]) << "*" << z[i][j] 
 
969
                            // << " + " << hat_node[idof][0] << "*" << z[i+1][j] << endl;
 
970
                        z_at_dof=(1-hat_node[idof][0])* z[i+1] + hat_node[idof][0] *z[i];
 
971
                        devy_at_dof=sqr(1-hat_node[idof][0])*z[i+1]/2*h[i] 
 
972
                                   -sqr(hat_node[idof][0])*z[i]/2*h[i] 
 
973
                                   +w[i];
 
974
                        n1=point(devy_at_dof[1],-devy_at_dof[0]);
 
975
                     }
 
976
                   kh.at(dofs[idof])= -1./(sqr(sqr(devy_at_dof[0])+sqr(devy_at_dof[1])))
 
977
                                *(z_at_dof[1]*devy_at_dof[0]-z_at_dof[0]*devy_at_dof[1])
 
978
                                *n1[j];
 
979
                  
 
980
                 }
 
981
  }
 
982
  return kh;
 
983
}
 
984
 
 
985
field
 
986
georep::plane_curvature_quadratic (const space& Nh, const domain& d, 
 
987
        const interface& bc) const
 
988
{
 
989
    check_macro(Nh.get_geo().dimension()==2,"Plane curvature is for 2D geometries only");
 
990
    check_macro(Nh.n_component()==2,"2-component vector space needed");
 
991
    // Sort contiguous elements of the interface
 
992
    //  first, group them in contiguous stubs listing element #s that are merged
 
993
    //  together as the jigsaw progresses...
 
994
    sort_interface(d,bc);
 
995
 
 
996
    // get coordinates of Lagrange interpolation points on edge elements:
 
997
    vector<point> hat_node;
 
998
    domain::const_iterator bd=d.begin();
 
999
    Nh.get_basis().hat_node(*bd, hat_node);
 
1000
 
 
1001
    Vector<size_t>::const_iterator iAB=bc.begin();
 
1002
    Vector<size_t>::const_iterator iBC=bc.begin(); iBC++;
 
1003
    Vector<size_t>::const_iterator the_end=bc.end();
 
1004
    // Calculate interface curvature using the mean of the curvature of
 
1005
    // circles passing through points e--,e-,e+ (curv1) and e-,e+,e++ (curv2)
 
1006
    field kh(Nh, Float(0));
 
1007
    tiny_vector<size_type> dofs;
 
1008
    point vcurv1, vcurv2;
 
1009
    point AB0=_x.begin()[bd[*iAB][0]];
 
1010
    point AB1=_x.begin()[bd[*iAB][1]];
 
1011
    point BC0=_x.begin()[bd[*iBC][0]];
 
1012
    point BC1=_x.begin()[bd[*iBC][1]];
 
1013
    point AB, BC;
 
1014
    bool rev_AB;
 
1015
    // check that nodes are in the same order as elements, otherwise reverse
 
1016
    if (AB1==BC0) { rev_AB=false; AB=AB1-AB0; BC=BC1-BC0; }
 
1017
    else
 
1018
      if (AB0==BC0) { rev_AB=true; AB=AB0-AB1; BC=BC1-BC0; }
 
1019
      else
 
1020
        if (AB1==BC1) { rev_AB=false; AB=AB1-AB0; BC=BC0-BC1;  }
 
1021
        else
 
1022
         {
 
1023
            check_macro(AB0==BC1, "Fatal: non consecutive edges!");
 
1024
            rev_AB=true; AB=AB0-AB1; BC=BC0-BC1; 
 
1025
         }
 
1026
    if (bc.first_tangent_imposed())
 
1027
    {
 
1028
        vcurv1=vcurvature2d(AB-2*dot(AB,bc.first_normal())*bc.first_normal(), AB);
 
1029
    }
 
1030
    else
 
1031
        vcurv1=vcurvature2d(AB,BC);
 
1032
    for ( ; iBC!=the_end; iAB++, iBC++) 
 
1033
    {
 
1034
        BC0=_x.begin()[bd[*iBC][0]];
 
1035
        BC1=_x.begin()[bd[*iBC][1]];
 
1036
        // check that nodes are in the same order as elements, otherwise reverse
 
1037
        if (AB1==BC0) { rev_AB=false; AB=AB1-AB0; BC=BC1-BC0; }
 
1038
        else
 
1039
          if (AB0==BC0) { rev_AB=true; AB=AB0-AB1; BC=BC1-BC0; }
 
1040
          else
 
1041
            if (AB1==BC1) { rev_AB=false; AB=AB1-AB0; BC=BC0-BC1;  }
 
1042
            else
 
1043
             {
 
1044
                check_macro(AB0==BC1, "Fatal: non consecutive edges!");
 
1045
                rev_AB=true; AB=AB0-AB1; BC=BC0-BC1; 
 
1046
             }
 
1047
        vcurv2=vcurvature2d(AB,BC);
 
1048
        /*
 
1049
        check_macro(abs(norm(vcurv2)-abs(curvature2d(AB,BC)))<1e-12,
 
1050
                "curvature error " << AB << " U " << BC << " " << norm(vcurv2) 
 
1051
                << "-" << curvature2d(AB,BC) << "=" << norm(vcurv2)-abs(curvature2d(AB,BC)));
 
1052
        */
 
1053
        for (size_type i=0; i<Nh.n_component(); i++) 
 
1054
        {
 
1055
            // get nh dofs corresponding
 
1056
            Nh.set_dof((Nh.get_geo().begin())[*iAB],dofs,i);
 
1057
            // set dofs
 
1058
            for (size_type idof=0; idof<dofs.size(); idof++)
 
1059
            if (!rev_AB)
 
1060
              kh.at(dofs[idof])=(1-hat_node[idof][0])*vcurv1[i] + hat_node[idof][0]*vcurv2[i];
 
1061
            else
 
1062
              kh.at(dofs[idof])=(1-hat_node[idof][0])*vcurv2[i] + hat_node[idof][0]*vcurv1[i];
 
1063
        }
 
1064
        vcurv1=vcurv2;
 
1065
        AB0=BC0; AB1=BC1; AB=BC;
 
1066
    }
 
1067
    // don't need to check that nodes are in the same order as elements, as BC was
 
1068
    rev_AB=(AB==AB0-AB1);
 
1069
    if (bc.last_tangent_imposed())
 
1070
    {
 
1071
        vcurv2=vcurvature2d(AB, AB-2*dot(AB,bc.last_normal())*bc.last_normal());
 
1072
    }
 
1073
    else vcurv2=vcurv1;
 
1074
    for (size_type i=0; i<Nh.n_component(); i++)
 
1075
    {
 
1076
        // get nh dofs corresponding
 
1077
        Nh.set_dof((Nh.get_geo().begin())[*iAB],dofs,i);
 
1078
        // set dofs
 
1079
        for (size_type idof=0; idof<dofs.size(); idof++)
 
1080
        if (!rev_AB)
 
1081
          kh.at(dofs[idof])=(1-hat_node[idof][0])*vcurv1[i] + hat_node[idof][0]*vcurv2[i];
 
1082
        else
 
1083
          kh.at(dofs[idof])=(1-hat_node[idof][0])*vcurv2[i] + hat_node[idof][0]*vcurv1[i];
 
1084
    }
 
1085
    return kh;
 
1086
}
 
1087
 
 
1088
 
 
1089
field
 
1090
georep::axi_curvature (const space& Nh, const domain& d, 
 
1091
        const interface& bc) const
 
1092
{
 
1093
    check_macro(Nh.coordinate_system()=="rz"||Nh.coordinate_system()=="zr",
 
1094
        "Axisymmetric curvature is for 2D-axisymmetric geometries only");
 
1095
    check_macro(Nh.n_component()==2,"2-component vector space needed");
 
1096
    // Sort contiguous elements of the interface
 
1097
    //  first, group them in contiguous stubs listing element #s that are merged
 
1098
    //  together as the jigsaw progresses...
 
1099
    sort_interface(d,bc);
 
1100
 
 
1101
    point er;
 
1102
    if (Nh.coordinate_system()=="rz") { er=point(1,0);  }
 
1103
        else { er=point(0,1);  }
 
1104
    // fill in with zeros in case only part of domain d is bounding region
 
1105
    field nh(Nh, Float(0));
 
1106
 
 
1107
    // get coordinates of Lagrange interpolation points on edge elements:
 
1108
    vector<point> hat_node;
 
1109
    domain::const_iterator bd=d.begin();
 
1110
    Nh.get_basis().hat_node(*bd, hat_node);
 
1111
 
 
1112
    Vector<size_t>::const_iterator iAB=bc.begin();
 
1113
    Vector<size_t>::const_iterator iBC=bc.begin(); iBC++;
 
1114
    Vector<size_t>::const_iterator the_end=bc.end();
 
1115
    field kh(Nh, Float(0));
 
1116
    tiny_vector<size_type> dofs;
 
1117
    Float curv0, curv1, curv2;
 
1118
    point nA, nB;
 
1119
    point AB0=_x.begin()[bd[*iAB][0]];
 
1120
    point AB1=_x.begin()[bd[*iAB][1]];
 
1121
    point BC0=_x.begin()[bd[*iBC][0]];
 
1122
    point BC1=_x.begin()[bd[*iBC][1]];
 
1123
    point AB, BC;
 
1124
    bool rev_AB=false;
 
1125
    // check that nodes are in the same order as elements, otherwise reverse
 
1126
    if (AB1==BC0) { rev_AB=false; AB=AB1-AB0; BC=BC1-BC0; }
 
1127
    else
 
1128
      if (AB0==BC0) { rev_AB=true; AB=AB0-AB1; BC=BC1-BC0; }
 
1129
      else
 
1130
        if (AB1==BC1) { rev_AB=false; AB=AB1-AB0; BC=BC0-BC1;  }
 
1131
        else
 
1132
         {
 
1133
            check_macro(AB0==BC1, "Fatal: non consecutive edges!");
 
1134
            rev_AB=true; AB=AB0-AB1; BC=BC0-BC1; 
 
1135
         }
 
1136
    curv0=curv1=curvature2d_axi((AB0+AB1)/2, AB, er);
 
1137
    if (bc.first_tangent_imposed())
 
1138
        nA=normal2d(AB-2*dot(AB,bc.first_normal())*bc.first_normal(), AB);
 
1139
    else
 
1140
        nA=normal2d(AB,BC);
 
1141
    for ( ; iBC!=the_end; iAB++, iBC++) 
 
1142
    {
 
1143
        BC0=_x.begin()[bd[*iBC][0]];
 
1144
        BC1=_x.begin()[bd[*iBC][1]];
 
1145
        // check that nodes are in the same order as elements, otherwise reverse
 
1146
        if (AB1==BC0) { rev_AB=false; AB=AB1-AB0; BC=BC1-BC0; }
 
1147
        else
 
1148
          if (AB0==BC0) { rev_AB=true; AB=AB0-AB1; BC=BC1-BC0; }
 
1149
          else
 
1150
            if (AB1==BC1) { rev_AB=false; AB=AB1-AB0; BC=BC0-BC1;  }
 
1151
            else
 
1152
             {
 
1153
                check_macro(AB0==BC1, "Fatal: non consecutive edges!");
 
1154
                rev_AB=true; AB=AB0-AB1; BC=BC0-BC1; 
 
1155
             }
 
1156
        nB=normal2d(AB, BC);
 
1157
        curv2=curvature2d_axi((BC0+BC1)/2, BC, er);
 
1158
        for (size_type i=0; i<Nh.n_component(); i++) 
 
1159
        {
 
1160
            // get nh dofs corresponding
 
1161
            Nh.set_dof((Nh.get_geo().begin())[*iAB],dofs,i);
 
1162
            // set dofs
 
1163
            for (size_type idof=0; idof<dofs.size(); idof++)
 
1164
            if (!rev_AB)
 
1165
              kh.at(dofs[idof])=.5*(1-hat_node[idof][0])*(curv0+curv1)*nA[i]
 
1166
                + .5*hat_node[idof][0]*(curv1+curv2)*nB[i];
 
1167
            else
 
1168
              kh.at(dofs[idof])=.5*hat_node[idof][0]*(curv0+curv1)*nA[i]
 
1169
                + .5*(1-hat_node[idof][0])*(curv1+curv2)*nB[i];
 
1170
        }
 
1171
        curv0=curv1;
 
1172
        curv1=curv2;
 
1173
        AB0=BC0; AB1=BC1; AB=BC;
 
1174
        nA=nB;
 
1175
    }
 
1176
    // don't need to check that nodes are in the same order as elements, as BC was
 
1177
    rev_AB=(AB==AB0-AB1);
 
1178
    curv2=curv1;
 
1179
    if (bc.last_tangent_imposed())
 
1180
        nB=normal2d(AB, AB-2*dot(AB,bc.last_normal())*bc.last_normal());
 
1181
    else 
 
1182
        nB=nA;
 
1183
    for (size_type i=0; i<Nh.n_component(); i++)
 
1184
    {
 
1185
        // get nh dofs corresponding
 
1186
        Nh.set_dof((Nh.get_geo().begin())[*iAB],dofs,i);
 
1187
        // set dofs
 
1188
        for (size_type idof=0; idof<dofs.size(); idof++)
 
1189
            if (!rev_AB)
 
1190
              kh.at(dofs[idof])=.5*(1-hat_node[idof][0])*(curv0+curv1)*nA[i]
 
1191
                + .5*hat_node[idof][0]*(curv1+curv2)*nB[i];
 
1192
            else
 
1193
              kh.at(dofs[idof])=.5*hat_node[idof][0]*(curv0+curv1)*nA[i]
 
1194
                + .5*(1-hat_node[idof][0])*(curv1+curv2)*nB[i];
 
1195
    }
 
1196
    return kh;
 
1197
}
 
1198
 
 
1199
void
 
1200
geo::interface_process (const domain& d, const interface& bc,
 
1201
        geometric_event& event) const
 
1202
 {
 
1203
    if (!bc.initialized())
 
1204
     {
 
1205
        space Nh(*this, d, "P0");
 
1206
        data().get_normal(Nh, d, false, bc);
 
1207
     }
 
1208
    data().interface_process(d, bc, event);
 
1209
 }
 
1210
void
 
1211
georep::interface_process (const domain& d, const interface& bc,
 
1212
        geometric_event& event) const
 
1213
 {
 
1214
    sort_interface(d, bc);
 
1215
    vector<point> hat_node;
 
1216
    domain::const_iterator bd=d.begin();
 
1217
 
 
1218
    Vector<size_t>::const_iterator iAB=bc.begin();
 
1219
    Vector<size_t>::const_iterator iBC=bc.begin(); iBC++;
 
1220
    Vector<size_t>::const_iterator the_end=bc.end();
 
1221
    
 
1222
    point AB0=_x.begin()[bd[*iAB][0]];
 
1223
    point AB1=_x.begin()[bd[*iAB][1]];
 
1224
    point BC0=_x.begin()[bd[*iBC][0]];
 
1225
    point BC1=_x.begin()[bd[*iBC][1]];
 
1226
    point A, B;
 
1227
    cerr << "!\tA=" << AB0;
 
1228
    for ( ; iBC!=the_end; iAB++, iBC++) 
 
1229
    {
 
1230
        BC0=_x.begin()[bd[*iBC][0]];
 
1231
        BC1=_x.begin()[bd[*iBC][1]];
 
1232
        // check that nodes are in the same order as elements, otherwise reverse
 
1233
        if (AB1==BC0) { A=AB0; B=AB1; }
 
1234
        else
 
1235
          if (AB0==BC0) { A=AB1; B=AB0; }
 
1236
          else
 
1237
            if (AB1==BC1) { A=AB0; B=AB1;  }
 
1238
            else
 
1239
             {
 
1240
                check_macro(AB0==BC1, "Fatal: non consecutive edges!");
 
1241
                A=AB1; B=AB0; 
 
1242
             }
 
1243
        event(A,B); 
 
1244
        AB0=BC0; AB1=BC1; A=B;
 
1245
    }
 
1246
    cerr << " B=" << BC1 << endl;
 
1247
cerr << "@";
 
1248
    if (BC0==B) event(A,BC1); else event(A,BC0);
 
1249
 }