Loading [MathJax]/extensions/tex2jax.js
Plugins
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
damage_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2019: Ruhr University Bochum
3 * Authors: Andreas Vogel
4 *
5 * This file is part of UG4.
6 *
7 * UG4 is free software: you can redistribute it and/or modify it under the
8 * terms of the GNU Lesser General Public License version 3 (as published by the
9 * Free Software Foundation) with the following additional attribution
10 * requirements (according to LGPL/GPL v3 §7):
11 *
12 * (1) The following notice must be displayed in the Appropriate Legal Notices
13 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 *
15 * (2) The following notice must be displayed at a prominent place in the
16 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 *
18 * (3) The following bibliography is recommended for citation and must be
19 * preserved in all covered files:
20 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 * parallel geometric multigrid solver on hierarchically distributed grids.
22 * Computing and visualization in science 16, 4 (2013), 151-164"
23 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 * flexible software system for simulating pde based models on high performance
25 * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU Lesser General Public License for more details.
31 */
32
33#ifndef __SMALL_STRAIN_MECHANICS__DAMAGE_IMPL_H_
34#define __SMALL_STRAIN_MECHANICS__DAMAGE_IMPL_H_
35
36#include "damage.h"
37#include "common/math/math_vector_matrix/math_vector_functions.h"
38
39namespace ug{
40namespace SmallStrainMechanics{
41
42
43template <int dim>
44void AveragePositions( MathVector<dim>& vCenter,
45 const std::vector<MathVector<dim> >& vCornerCoords)
46{
47 vCenter = vCornerCoords[0];
48 for(size_t j = 1; j < vCornerCoords.size(); ++j)
49 {
50 vCenter += vCornerCoords[j];
51 }
52 vCenter *= 1./(number)( vCornerCoords.size());
53}
54
55
57// Collect Surface Neighbors
59
60
61template <typename TDomain>
62void CollectSurfaceNeighbors(
63 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF, // some dummy function
64 typename grid_dim_traits<TDomain::dim>::element_type* elem,
65 std::vector< typename grid_dim_traits<TDomain::dim>::element_type * >& vNeighbors)
66{
67 PROFILE_BEGIN_GROUP(CollectSurfaceNeighbors, "Small Strain Mech");
68
69 static const int dim = TDomain::dim;
70 typedef typename TDomain::grid_type TGrid;
71 typedef typename grid_dim_traits<dim>::element_type TElem;
72 typedef typename grid_dim_traits<dim>::side_type TSide;
73 typedef typename contrained_dim_traits<dim>::contrained_side_type TContrainedSide;
74 typedef typename contrained_dim_traits<dim>::contraining_side_type TContrainingSide;
75
76 typedef typename TDomain::position_accessor_type TPositionAccessor;
77
78 // get domain
79 SmartPtr<TDomain> domain = spF->domain();
80 SmartPtr<typename TDomain::grid_type> spGrid = domain->grid();
81 typename TDomain::grid_type& grid = *spGrid;
82
83 // get dof distribution
84 SmartPtr<DoFDistribution> dd = spF->dd();
85
86 // get iterators
87 typename DoFDistribution::traits<TElem>::iterator iter, iterEnd;
88 iter = dd->begin<TElem>(SurfaceView::ALL); // SurfaceView::MG_ALL
89 iterEnd = dd->end<TElem>(SurfaceView::ALL);
90
91
92 // clear container
93 vNeighbors.clear();
94
95 // get all sides in order
96 typename TGrid::template traits<TSide>::secure_container vSide;
97 grid.associated_elements_sorted(vSide, elem);
98
99 // begin marking
100 grid.begin_marking();
101
102 // mark the elem
103 grid.mark(elem);
104
106 // loop all sides
108 for(size_t s = 0; s < vSide.size(); ++s){
109
110 // get all connected elements
111 typename TGrid::template traits<TElem>::secure_container vElemOfSide;
112 grid.associated_elements(vElemOfSide, vSide[s]);
113
114 // if no elem found: internal error
115 if(vElemOfSide.size() == 0)
116 UG_THROW("Huh, why no element? Should be at least the considered elem itself");
117
119 // handle sides with only one element
120 // ( if only one element for side, it must be a boundary element or contrained)
122 if(vElemOfSide.size() == 1){
123
125 // handle constraint edge
126 //
128 if(vSide[s]->is_constrained()){
129
130 TContrainedSide* cSide = dynamic_cast<TContrainedSide*>(vSide[s]);
131 TSide* constrainingSide = dynamic_cast<TSide*>(cSide->get_constraining_object());
132
133 typename TGrid::template traits<TElem>::secure_container vElemOfContrainingSide;
134 grid.associated_elements(vElemOfContrainingSide, constrainingSide);
135 if(vElemOfContrainingSide.size() != 2) UG_THROW("Huh, should be "<<2<<" at constraining side");
136
137 for(size_t nbr = 0; nbr < vElemOfContrainingSide.size(); ++nbr){
138
139 TElem* neighborElem = vElemOfContrainingSide[nbr];
140 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
141
142 grid.mark(neighborElem);
143 vNeighbors.push_back(neighborElem);
144 }
145 }
146 }
148 // handle hanging sides with 2 elements
149 // ( coupled to all fine neighbor elems )
151 else if(vSide[s]->is_constraining()){
152
153 TContrainingSide* cSide = dynamic_cast<TContrainingSide*>(vSide[s]);
154 const size_t numConstrained = cSide->template num_constrained<TSide>();
155
156 for(size_t cs = 0; cs < numConstrained; ++cs){
157
158 TSide* constrainedSide = cSide->template constrained<TSide>(cs);
159
160 // neighbor elem
161 typename TGrid::template traits<TElem>::secure_container vElemOfContrainedSide;
162 grid.associated_elements(vElemOfContrainedSide, constrainedSide);
163 if(vElemOfContrainedSide.size() != 1) UG_THROW("Huh, should be 1 at constrained side");
164 TElem* neighborElem = vElemOfContrainedSide[0];
165 grid.mark(neighborElem);
166 vNeighbors.push_back(neighborElem);
167
168 }
169
170 }
172 // regular refined sides
173 // (find all direct face neighbors)
175 else{
176 // if more than 2 elem found: internal error
177 if(vElemOfSide.size() != 2)
178 UG_THROW("Huh, why more than 2 elements of side?");
179
180 for(size_t eos = 0; eos < vElemOfSide.size(); ++eos){
181
182 // neighbor elem
183 TElem* neighborElem = vElemOfSide[eos];
184
185 // check
186 if(grid.template num_children<TElem,TElem>(neighborElem) > 0)
187 UG_THROW("Huh, why not on top level?");
188
189 // skip self
190 if(neighborElem == elem) continue;
191 grid.mark(neighborElem);
192 vNeighbors.push_back(neighborElem);
193 }
194 }
195
196 } // end loop sides
197
199 // loop all sides (FOR DIAG COUPLING)
201 for(size_t s = 0; s < vSide.size(); ++s){
202
203 // get all connected elements
204 typename TGrid::template traits<TElem>::secure_container vElemOfSide;
205 grid.associated_elements(vElemOfSide, vSide[s]);
206
207 // if no elem found: internal error
208 if(vElemOfSide.size() == 0)
209 UG_THROW("Huh, why no element? Should be at least the considered elem itself");
210
211
213 // diagonal elements (only coupled via a common vertex)
214 // (find all vertex coupled neighbors)
216
217 // loop vertices of side
218 Vertex* const* vVertex = vSide[s]->vertices();
219 const size_t numVertex = vSide[s]->num_vertices();
220
221
222 for(size_t vrt = 0; vrt < numVertex; ++vrt)
223 {
224 // collect all diag leaf-elems on same level
225 typename TGrid::template traits<TElem>::secure_container vVertexNeighbor;
226 grid.associated_elements(vVertexNeighbor, vVertex[vrt]);
227 for(size_t eov = 0; eov < vVertexNeighbor.size(); ++eov)
228 {
229 TElem* neighborElem = vVertexNeighbor[eov];
230 if(grid.is_marked(neighborElem)) continue;
231 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
232
233 grid.mark(neighborElem);
234 vNeighbors.push_back(neighborElem);
235 }
236
237 // collect all diag leaf-elems on finer level
238 if(grid.template num_children<Vertex,Vertex>(vVertex[vrt]) > 0){
239
240 if(grid.template num_children<Vertex,Vertex>(vVertex[vrt]) != 1)
241 UG_THROW("Huh, why more than one vertex child?")
242
243 Vertex* fineVrt = grid.template get_child<Vertex,Vertex>(vVertex[vrt], 0);
244
245 typename TGrid::template traits<TElem>::secure_container vVertexNeighbor;
246 grid.associated_elements(vVertexNeighbor, fineVrt);
247 for(size_t eov = 0; eov < vVertexNeighbor.size(); ++eov)
248 {
249 TElem* neighborElem = vVertexNeighbor[eov];
250 if(grid.is_marked(neighborElem)) continue;
251 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
252
253 grid.mark(neighborElem);
254 vNeighbors.push_back(neighborElem);
255 }
256
257 }
258
259 // collect all diag leaf-elems on coarser level
260 if(grid.get_parent(vVertex[vrt]) != 0){
261
262
263 Vertex* coarseVrt = dynamic_cast<Vertex*>(grid.get_parent(vVertex[vrt]));
264
265 if(coarseVrt != 0){
266
267 typename TGrid::template traits<TElem>::secure_container vVertexNeighbor;
268 grid.associated_elements(vVertexNeighbor, coarseVrt);
269 for(size_t eov = 0; eov < vVertexNeighbor.size(); ++eov)
270 {
271 TElem* neighborElem = vVertexNeighbor[eov];
272 if(grid.is_marked(neighborElem)) continue;
273 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
274
275 grid.mark(neighborElem);
276 vNeighbors.push_back(neighborElem);
277 }
278 }
279 }
280
281
282 }
283 } // end side loop
284
285
286 // end marking
287 grid.end_marking();
288
289}
290
292// CollectStencilNeighbors_NeumannZeroBND_IndexAndDistance
294
295template <typename TDomain>
296void CollectStencilNeighbors_NeumannZeroBND_IndexAndDistance
297(
298 std::vector< typename grid_dim_traits<TDomain::dim>::element_type* >& vElem,
299 std::vector<size_t>& vIndex,
300 std::vector< MathVector<TDomain::dim> >& vDistance,
301 typename grid_dim_traits<TDomain::dim>::element_type* elem,
302 typename TDomain::grid_type& grid,
303 typename TDomain::position_accessor_type& aaPos,
304 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF, bool fillElemSizeIntoVector
305)
306{
307 PROFILE_BEGIN_GROUP(DamageFunctionUpdater_CollectStencilNeighbors, "Small Strain Mech");
308
309 static const int dim = TDomain::dim;
310 typedef typename TDomain::grid_type TGrid;
311 typedef typename grid_dim_traits<dim>::element_type TElem;
312 typedef typename grid_dim_traits<dim>::side_type TSide;
313 typedef typename contrained_dim_traits<dim>::contrained_side_type TContrainedSide;
314 typedef typename contrained_dim_traits<dim>::contraining_side_type TContrainingSide;
315 typedef typename TDomain::position_accessor_type TPositionAccessor;
316
317
318 //static const int numNeighborsToFind = 2*dim + (dim * (dim-1)) / 2;
319 const size_t fct = 0;
320
321 vElem.clear();
322 vIndex.clear();
323 vDistance.clear();
324
325 // get vertices of element
326 Vertex* const* vVertex = elem->vertices();
327 const size_t numVertex = elem->num_vertices();
328
329 // corner coordinates
330 std::vector<MathVector<dim> > vCornerCoords;
331 vCornerCoords.resize(numVertex);
332 for(size_t vrt = 0; vrt < numVertex; ++vrt){
333 vCornerCoords[vrt] = aaPos[ vVertex[vrt] ];
334 }
335
336 // element midpoint
337 MathVector<dim> ElemMidPoint;
338 AveragePositions<dim>(ElemMidPoint, vCornerCoords);
339
340// UG_LOG("############## Element with midpoint " << ElemMidPoint << " ################ \n")
341
342 if(fillElemSizeIntoVector){
343
344 // add index
345 std::vector<DoFIndex> ind;
346 if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
347 const size_t i = ind[0][0];
348
349 // element volume
350 (*spF)[i] = ElementSize<dim>(elem->reference_object_id(), &vCornerCoords[0]);
351 }
352
353
354 // get all sides in order
355 typename TGrid::template traits<TSide>::secure_container vSide;
356 grid.associated_elements_sorted(vSide, elem);
357
358 // begin marking
359 grid.begin_marking();
360
361 // mark the elem
362 grid.mark(elem);
363
365 // loop all sides
367 for(size_t s = 0; s < vSide.size(); ++s){
368
369 // get all connected elements
370 typename TGrid::template traits<TElem>::secure_container vElemOfSide;
371 grid.associated_elements(vElemOfSide, vSide[s]);
372
373 // if no elem found: internal error
374 if(vElemOfSide.size() == 0)
375 UG_THROW("Huh, why no element?");
376
378 // handle sides with only one element
379 // ( if only one element for side, it must be a boundary element or contrained)
381 if(vElemOfSide.size() == 1){
382
384 // handle constraint edge
386 if(vSide[s]->is_constrained()){
387
388 TContrainedSide* cSide = dynamic_cast<TContrainedSide*>(vSide[s]);
389 TSide* constrainingSide = dynamic_cast<TSide*>(cSide->get_constraining_object());
390
391 typename TGrid::template traits<TElem>::secure_container vElemOfContrainingSide;
392 grid.associated_elements(vElemOfContrainingSide, constrainingSide);
393 if(vElemOfContrainingSide.size() != 2) UG_THROW("Huh, should be "<<2<<" at constraining side");
394
395 for(size_t nbr = 0; nbr < vElemOfContrainingSide.size(); ++nbr){
396
397 TElem* neighborElem = vElemOfContrainingSide[nbr];
398 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
399
400 grid.mark(neighborElem);
401 vElem.push_back(neighborElem);
402
403 // add index
404 std::vector<DoFIndex> ind;
405 if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs");
406 vIndex.push_back(ind[0][0]);
407
408 // add distance
409 vDistance.resize(vDistance.size()+1);
410 CollectCornerCoordinates(vCornerCoords, *neighborElem, aaPos);
411 AveragePositions<dim>(vDistance.back(), vCornerCoords);
412 VecScaleAdd(vDistance.back(), 1.0, ElemMidPoint, -1.0, vDistance.back());
413 }
414
415 }
416
418 // handle mirror elements
419 // ( if only one element for side, it must be a boundary element)
421 else {
422 vElem.push_back(vElemOfSide[0]);
423
424 // add index
425 std::vector<DoFIndex> ind;
426 if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
427 vIndex.push_back(ind[0][0]);
428
429 // add distance
430 MathVector<dim> n;
431 CollectCornerCoordinates(vCornerCoords, *vSide[s], aaPos);
432 ElementNormal<dim>(vSide[s]->reference_object_id(), n, &vCornerCoords[0]);
433
434 vDistance.resize(vDistance.size()+1);
435 ProjectPointToPlane(vDistance.back(), ElemMidPoint, vCornerCoords[0], n);
436 VecScaleAdd(vDistance.back(), 2.0, ElemMidPoint, -2.0, vDistance.back());
437// UG_LOG("is boundary, with vDistance: " << vDistance.back() << "\n");
438 }
439 }
441 // handle hanging sides with 2 elements
442 // ( coupled to all fine neighbor elems )
444 else if(vSide[s]->is_constraining()){
445
446 TContrainingSide* cSide = dynamic_cast<TContrainingSide*>(vSide[s]);
447 const size_t numConstrained = cSide->template num_constrained<TSide>();
448
449 for(size_t cs = 0; cs < numConstrained; ++cs){
450
451 TSide* constrainedSide = cSide->template constrained<TSide>(cs);
452
453 // neighbor elem
454 typename TGrid::template traits<TElem>::secure_container vElemOfContrainedSide;
455 grid.associated_elements(vElemOfContrainedSide, constrainedSide);
456 if(vElemOfContrainedSide.size() != 1) UG_THROW("Huh, should be 1 at constrained side");
457 TElem* neighborElem = vElemOfContrainedSide[0];
458 grid.mark(neighborElem);
459 vElem.push_back(neighborElem);
460
461 // add index
462 std::vector<DoFIndex> ind;
463 if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs");
464 vIndex.push_back(ind[0][0]);
465
466 // add distance
467 vDistance.resize(vDistance.size()+1);
468 CollectCornerCoordinates(vCornerCoords, *neighborElem, aaPos);
469 AveragePositions<dim>(vDistance.back(), vCornerCoords);
470 VecScaleAdd(vDistance.back(), 1.0, ElemMidPoint, -1.0, vDistance.back());
471 }
472
473
474 }
476 // regular refined sides
477 // (find all direct face neighbors)
479 else{
480 // if more than 2 elem found: internal error
481 if(vElemOfSide.size() != 2)
482 UG_THROW("Huh, why more than 2 elements of side?");
483
484 for(size_t eos = 0; eos < vElemOfSide.size(); ++eos){
485
486 // neighbor elem
487 TElem* neighborElem = vElemOfSide[eos];
488
489 // skip self
490 if(neighborElem == elem) continue;
491
492 // add neighbor
493 grid.mark(neighborElem);
494 vElem.push_back(neighborElem);
495
496 // add index
497 std::vector<DoFIndex> ind;
498 if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs");
499 vIndex.push_back(ind[0][0]);
500
501 // add distance
502 vDistance.resize(vDistance.size()+1);
503 CollectCornerCoordinates(vCornerCoords, *neighborElem, aaPos);
504 AveragePositions<dim>(vDistance.back(), vCornerCoords);
505 VecScaleAdd(vDistance.back(), 1.0, ElemMidPoint, -1.0, vDistance.back());
506// UG_LOG("has neighbor, c" << vDistance.back() << "\n");
507 }
508 }
509 }
510
512 // add additional elements (via vertex coupling)
514 for(size_t vrt = 0; vrt < numVertex; ++vrt)
515 {
516 typename TGrid::template traits<TElem>::secure_container vVertexNeighbor;
517 grid.associated_elements(vVertexNeighbor, vVertex[vrt]);
518
519 // collect all diag leaf-elems on same level
520 for(size_t eov = 0; eov < vVertexNeighbor.size(); ++eov)
521 {
522 TElem* neighborElem = vVertexNeighbor[eov];
523 if(grid.is_marked(neighborElem)) continue;
524 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
525
526 grid.mark(neighborElem);
527 vElem.push_back(neighborElem);
528
529 // add index
530 std::vector<DoFIndex> ind;
531 if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs");
532 vIndex.push_back(ind[0][0]);
533
534 // add distance
535 vDistance.resize(vDistance.size()+1);
536 CollectCornerCoordinates(vCornerCoords, *neighborElem, aaPos);
537 AveragePositions<dim>(vDistance.back(), vCornerCoords);
538 VecScaleAdd(vDistance.back(), 1.0, ElemMidPoint, -1.0, vDistance.back());
539 }
540
541
542 // collect all diag leaf-elems on finer level
543 if(grid.template num_children<Vertex,Vertex>(vVertex[vrt]) > 0){
544
545 if(grid.template num_children<Vertex,Vertex>(vVertex[vrt]) != 1)
546 UG_THROW("Huh, why more than one vertex child?")
547
548 Vertex* fineVrt = grid.template get_child<Vertex,Vertex>(vVertex[vrt], 0);
549
550 typename TGrid::template traits<TElem>::secure_container vVertexNeighbor;
551 grid.associated_elements(vVertexNeighbor, fineVrt);
552 for(size_t eov = 0; eov < vVertexNeighbor.size(); ++eov)
553 {
554 TElem* neighborElem = vVertexNeighbor[eov];
555 if(grid.is_marked(neighborElem)) continue;
556 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
557
558 grid.mark(neighborElem);
559 vElem.push_back(neighborElem);
560
561 // add index
562 std::vector<DoFIndex> ind;
563 if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs");
564 vIndex.push_back(ind[0][0]);
565
566 // add distance
567 vDistance.resize(vDistance.size()+1);
568 CollectCornerCoordinates(vCornerCoords, *neighborElem, aaPos);
569 AveragePositions<dim>(vDistance.back(), vCornerCoords);
570 VecScaleAdd(vDistance.back(), 1.0, ElemMidPoint, -1.0, vDistance.back());
571 }
572
573 }
574
575 // collect all diag leaf-elems on coarser level
576 if(grid.get_parent(vVertex[vrt]) != 0){
577
578 Vertex* coarseVrt = dynamic_cast<Vertex*>(grid.get_parent(vVertex[vrt]));
579
580 if(coarseVrt != 0){
581
582 typename TGrid::template traits<TElem>::secure_container vVertexNeighbor;
583 grid.associated_elements(vVertexNeighbor, coarseVrt);
584 for(size_t eov = 0; eov < vVertexNeighbor.size(); ++eov)
585 {
586 TElem* neighborElem = vVertexNeighbor[eov];
587 if(grid.is_marked(neighborElem)) continue;
588 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
589
590 grid.mark(neighborElem);
591 vElem.push_back(neighborElem);
592
593 // add index
594 std::vector<DoFIndex> ind;
595 if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs");
596 vIndex.push_back(ind[0][0]);
597
598 // add distance
599 vDistance.resize(vDistance.size()+1);
600 CollectCornerCoordinates(vCornerCoords, *neighborElem, aaPos);
601 AveragePositions<dim>(vDistance.back(), vCornerCoords);
602 VecScaleAdd(vDistance.back(), 1.0, ElemMidPoint, -1.0, vDistance.back());
603 }
604 }
605 }
606
607
608 }
609
610 // end marking
611 grid.end_marking();
612}
613
614
616// By partial integration
618
619
620template <typename TDomain>
621void InitLaplacian_PartialIntegration(
622 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF, // some dummy function
623 std::vector< std::vector< number > >& vStencil,
624 std::vector< std::vector<size_t> >& vIndex,
625 int quadRuleType, bool fillElemSizeIntoVector)
626{
627 PROFILE_BEGIN_GROUP(InitLaplacian_PartialIntegration, "Small Strain Mech");
628
629 static const int dim = TDomain::dim;
630 typedef typename TDomain::grid_type TGrid;
631 typedef typename grid_dim_traits<dim>::element_type TElem;
632 typedef typename grid_dim_traits<dim>::side_type TSide;
633 typedef typename contrained_dim_traits<dim>::contrained_side_type TContrainedSide;
634 typedef typename contrained_dim_traits<dim>::contraining_side_type TContrainingSide;
635
636 typedef typename TDomain::position_accessor_type TPositionAccessor;
637
638 if(dim == 3)
639 UG_THROW("This implementation is 2d only, currently. Handle vertex neighbors properly in 3d...");
640
641 const size_t fct = 0;
642
643 number sideWeight;
644 number vertexWeight;
645
646 if(quadRuleType == 1) // Midpoint
647 {
648 sideWeight = 1.0;
649 vertexWeight = 0.0;
650 }
651 else if (quadRuleType == 2) // Simpson's
652 {
653 sideWeight = 4.0/6.0;
654 vertexWeight = 1.0/6.0;
655 }
656 else
657 UG_THROW("InitLaplacian_PartialIntegration: wrong quad rule")
658
659
660 // get domain
661 SmartPtr<TDomain> domain = spF->domain();
662 SmartPtr<typename TDomain::grid_type> spGrid = domain->grid();
663 typename TDomain::grid_type& grid = *spGrid;
664 typename TDomain::position_accessor_type& aaPos = domain->position_accessor();
665
666 // get dof distribution
667 SmartPtr<DoFDistribution> dd = spF->dd();
668
669 // get iterators
670 typename DoFDistribution::traits<TElem>::iterator iter, iterEnd;
671 iter = dd->begin<TElem>(SurfaceView::ALL); // SurfaceView::MG_ALL
672 iterEnd = dd->end<TElem>(SurfaceView::ALL);
673
674 // resize result vectors
675 const size_t numIndex = spF->num_dofs();
676 vIndex.resize(numIndex);
677 vStencil.resize(numIndex);
678
679 // storage (for reuse)
680 std::vector<MathVector<dim> > vCornerCoords, vNbrCornerCoords;
681 MathVector<dim> ElemMidPoint, Normal, Distance;
682
683 // loop all vertices
684 for(;iter != iterEnd; ++iter)
685 {
686 // get vertex
687 TElem* elem = *iter;
688
689 // store index
690 std::vector<DoFIndex> ind;
691 if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
692 const size_t i = ind[0][0];
693
694 // add self-coupling
695 vIndex[i].clear();
696 vIndex[i].push_back(i);
697
698 // reset stencil
699 vStencil[i].clear();
700 vStencil[i].resize(1, 0.0);
701
702 // get vertices of element
703 Vertex* const* vVertex = elem->vertices();
704 const size_t numVertex = elem->num_vertices();
705
706 // corner coordinates
707 vCornerCoords.resize(numVertex);
708 for(size_t vrt = 0; vrt < numVertex; ++vrt){
709 vCornerCoords[vrt] = aaPos[ vVertex[vrt] ];
710 }
711
712 // element volume
713 const number vol = ElementSize<dim>(elem->reference_object_id(), &vCornerCoords[0]);
714
715 if(fillElemSizeIntoVector)
716 (*spF)[i] = vol;
717
718 // element midpoint
719 AveragePositions<dim>(ElemMidPoint, vCornerCoords);
720
721 // get all sides in order
722 typename TGrid::template traits<TSide>::secure_container vSide;
723 grid.associated_elements_sorted(vSide, elem);
724
725 // begin marking
726 grid.begin_marking();
727
728 // mark the elem
729 grid.mark(elem);
730
732 // loop all sides
734 for(size_t s = 0; s < vSide.size(); ++s){
735
736 // side normal
737 SideNormal<dim>(elem->reference_object_id(), Normal, s, &vCornerCoords[0]);
738
739 // get all connected elements
740 typename TGrid::template traits<TElem>::secure_container vElemOfSide;
741 grid.associated_elements(vElemOfSide, vSide[s]);
742
743 // if no elem found: internal error
744 if(vElemOfSide.size() == 0)
745 UG_THROW("Huh, why no element? Should be at least the considered elem itself");
746
748 // handle sides with only one element
749 // ( if only one element for side, it must be a boundary element or contrained)
751 if(vElemOfSide.size() == 1){
752
754 // handle constraint edge
755 //
757 if(vSide[s]->is_constrained()){
758
759 TContrainedSide* cSide = dynamic_cast<TContrainedSide*>(vSide[s]);
760 TSide* constrainingSide = dynamic_cast<TSide*>(cSide->get_constraining_object());
761
762 typename TGrid::template traits<TElem>::secure_container vElemOfContrainingSide;
763 grid.associated_elements(vElemOfContrainingSide, constrainingSide);
764 if(vElemOfContrainingSide.size() != 2) UG_THROW("Huh, should be 2 at constraining side");
765
766 for(size_t nbr = 0; nbr < vElemOfContrainingSide.size(); ++nbr){
767
768 TElem* neighborElem = vElemOfContrainingSide[nbr];
769 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
770
771 grid.mark(neighborElem);
772
773 // add distance
774 CollectCornerCoordinates(vNbrCornerCoords, *neighborElem, aaPos);
775 AveragePositions<dim>(Distance, vNbrCornerCoords);
776 VecScaleAdd(Distance, 1.0, ElemMidPoint, -1.0, Distance);
777
778 const number cpl = (sideWeight) * VecDot(Normal, Distance) / ( VecTwoNormSq(Distance) * vol);
779 vStencil[i].push_back( -cpl );
780 vStencil[i][0] += cpl;
781
782 std::vector<DoFIndex> ind;
783 if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs");
784 vIndex[i].push_back(ind[0][0]);
785
786 }
787 }
789 // handle mirror elements
790 // ( if only one element for side, it must be a boundary element)
792 else {
793
794 ProjectPointToPlane(Distance, ElemMidPoint, aaPos[ (vSide[s]->vertex(0)) ], Normal);
795 VecScaleAdd(Distance, 2.0, ElemMidPoint, -2.0, Distance);
796
797 const number cpl = VecDot(Normal, Distance) / ( VecTwoNormSq(Distance) * vol);
798 vStencil[i].push_back( -cpl );
799 vStencil[i][0] += cpl;
800
801 std::vector<DoFIndex> ind;
802 if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
803 vIndex[i].push_back(ind[0][0]);
804 }
805 }
807 // handle hanging sides with 2 elements
808 // ( coupled to all fine neighbor elems )
810 else if(vSide[s]->is_constraining()){
811
812 TContrainingSide* cSide = dynamic_cast<TContrainingSide*>(vSide[s]);
813 const size_t numConstrained = cSide->template num_constrained<TSide>();
814
815 for(size_t cs = 0; cs < numConstrained; ++cs){
816
817 TSide* constrainedSide = cSide->template constrained<TSide>(cs);
818
819 // neighbor elem
820 typename TGrid::template traits<TElem>::secure_container vElemOfContrainedSide;
821 grid.associated_elements(vElemOfContrainedSide, constrainedSide);
822 if(vElemOfContrainedSide.size() != 1) UG_THROW("Huh, should be 1 at constrained side");
823 TElem* neighborElem = vElemOfContrainedSide[0];
824 grid.mark(neighborElem);
825
826 // add distance
827 CollectCornerCoordinates(vNbrCornerCoords, *neighborElem, aaPos);
828 AveragePositions<dim>(Distance, vNbrCornerCoords);
829 VecScaleAdd(Distance, 1.0, ElemMidPoint, -1.0, Distance);
830
831 const number cpl = (sideWeight / numConstrained) * VecDot(Normal, Distance) / ( VecTwoNormSq(Distance) * vol);
832 vStencil[i].push_back( -cpl );
833 vStencil[i][0] += cpl;
834
835 std::vector<DoFIndex> ind;
836 if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs");
837 vIndex[i].push_back(ind[0][0]);
838 }
839
840 }
842 // regular refined sides
843 // (find all direct face neighbors)
845 else{
846 // if more than 2 elem found: internal error
847 if(vElemOfSide.size() != 2)
848 UG_THROW("Huh, why more than 2 elements of side?");
849
850 for(size_t eos = 0; eos < vElemOfSide.size(); ++eos){
851
852 // neighbor elem
853 TElem* neighborElem = vElemOfSide[eos];
854
855 // check
856 if(grid.template num_children<TElem,TElem>(neighborElem) > 0)
857 UG_THROW("Huh, why not on top level?");
858
859 // skip self
860 if(neighborElem == elem) continue;
861 grid.mark(neighborElem);
862
863 // add distance
864 CollectCornerCoordinates(vNbrCornerCoords, *neighborElem, aaPos);
865 AveragePositions<dim>(Distance, vNbrCornerCoords);
866 VecScaleAdd(Distance, 1.0, ElemMidPoint, -1.0, Distance);
867
868 const number cpl = sideWeight * VecDot(Normal, Distance) / ( VecTwoNormSq(Distance) * vol);
869 vStencil[i].push_back( -cpl );
870 vStencil[i][0] += cpl;
871
872 std::vector<DoFIndex> ind;
873 if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs");
874 vIndex[i].push_back(ind[0][0]);
875 }
876 }
877
878 } // end loop sides
879
881 // loop all sides (FOR DIAG COUPLING)
883 if(vertexWeight != 0.0){
884
885 for(size_t s = 0; s < vSide.size(); ++s){
886
887 // side normal
888 SideNormal<dim>(elem->reference_object_id(), Normal, s, &vCornerCoords[0]);
889
890 // get all connected elements
891 typename TGrid::template traits<TElem>::secure_container vElemOfSide;
892 grid.associated_elements(vElemOfSide, vSide[s]);
893
894 // if no elem found: internal error
895 if(vElemOfSide.size() == 0)
896 UG_THROW("Huh, why no element? Should be at least the considered elem itself");
897
898
900 // diagonal elements (only coupled via a common vertex)
901 // (find all vertex coupled neighbors)
903
904 // loop vertices of side
905 Vertex* const* vVertex = vSide[s]->vertices();
906 const size_t numVertex = vSide[s]->num_vertices();
907
908
909 for(size_t vrt = 0; vrt < numVertex; ++vrt)
910 {
911 std::vector<TElem*> vDiagNeighbors;
912
913 // collect all diag leaf-elems on same level
914 typename TGrid::template traits<TElem>::secure_container vVertexNeighbor;
915 grid.associated_elements(vVertexNeighbor, vVertex[vrt]);
916 for(size_t eov = 0; eov < vVertexNeighbor.size(); ++eov)
917 {
918 TElem* neighborElem = vVertexNeighbor[eov];
919 if(grid.is_marked(neighborElem)) continue;
920 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
921
922 //UG_LOG("Over Same level\n")
923
924 //grid.mark(neighborElem);
925
926 vDiagNeighbors.push_back(neighborElem);
927 }
928
929 // collect all diag leaf-elems on finer level
930 if(grid.template num_children<Vertex,Vertex>(vVertex[vrt]) > 0){
931
932 if(grid.template num_children<Vertex,Vertex>(vVertex[vrt]) != 1)
933 UG_THROW("Huh, why more than one vertex child?")
934
935 Vertex* fineVrt = grid.template get_child<Vertex,Vertex>(vVertex[vrt], 0);
936
937 typename TGrid::template traits<TElem>::secure_container vVertexNeighbor;
938 grid.associated_elements(vVertexNeighbor, fineVrt);
939 for(size_t eov = 0; eov < vVertexNeighbor.size(); ++eov)
940 {
941 TElem* neighborElem = vVertexNeighbor[eov];
942 if(grid.is_marked(neighborElem)) continue;
943 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
944
945 //UG_LOG("Over Children\n")
946
947 //grid.mark(neighborElem);
948
949 vDiagNeighbors.push_back(neighborElem);
950 }
951
952 }
953
954 // collect all diag leaf-elems on coarser level
955 if(grid.get_parent(vVertex[vrt]) != 0){
956
957
958 Vertex* coarseVrt = dynamic_cast<Vertex*>(grid.get_parent(vVertex[vrt]));
959
960 if(coarseVrt != 0){
961
962 typename TGrid::template traits<TElem>::secure_container vVertexNeighbor;
963 grid.associated_elements(vVertexNeighbor, coarseVrt);
964 for(size_t eov = 0; eov < vVertexNeighbor.size(); ++eov)
965 {
966 TElem* neighborElem = vVertexNeighbor[eov];
967 if(grid.is_marked(neighborElem)) continue;
968 if(grid.template num_children<TElem,TElem>(neighborElem) > 0) continue;
969
970 //UG_LOG("Over Parent\n")
971
972 //grid.mark(neighborElem);
973
974 vDiagNeighbors.push_back(neighborElem);
975 }
976 }
977 }
978
979 // add contribution
980 const size_t numDiagElem = vDiagNeighbors.size();
981 for(size_t diag = 0; diag < numDiagElem; ++diag){
982
983 TElem* neighborElem = vDiagNeighbors[diag];
984
985 // add distance
986 CollectCornerCoordinates(vNbrCornerCoords, *neighborElem, aaPos);
987 AveragePositions<dim>(Distance, vNbrCornerCoords);
988 VecScaleAdd(Distance, 1.0, ElemMidPoint, -1.0, Distance);
989
990
991 //UG_LOG("Elem: "<<i<<", side: "<<s<<", vrt: "<<vrt<<", Distance: "<<Distance<<", Normal: "<<Normal<<"\n");
992
993 std::vector<DoFIndex> ind;
994 if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs");
995 const size_t j = ind[0][0];
996
997 size_t k = 0;
998 for(; k < vIndex[i].size(); ++k){
999 if(vIndex[i][k] == j)
1000 break;
1001 }
1002 if(k == vIndex[i].size())
1003 {
1004 vIndex[i].push_back(j);
1005 vStencil[i].push_back(0.0);
1006 // k == vIndex[i].size();
1007 }
1008
1009 const number cpl = (vertexWeight / numDiagElem) * VecDot(Normal, Distance) / ( VecTwoNormSq(Distance) * vol);
1010 vStencil[i][k] -= cpl;
1011 vStencil[i][0] += cpl;
1012 }
1013 }
1014 } // end side loop
1015
1016 } // end diag-elem block
1017
1018 // end marking
1019 grid.end_marking();
1020
1021 } // end element loop
1022
1023 static int call = 0; call++;
1024 //write_stencil_matrix_debug(spF, "Stencil", call);
1025
1026}
1027
1028
1030// By taylor expansion
1032
1033
1034template <typename TDomain>
1035void InitLaplacian_TaylorExpansion(
1036 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF,
1037 std::vector< std::vector< number > >& vStencil,
1038 std::vector< std::vector<size_t> >& vIndex, bool fillElemSizeIntoVector)
1039{
1040 PROFILE_BEGIN_GROUP(DamageFunctionUpdater_init_TaylorExpansion, "Small Strain Mech");
1041
1042 static const int dim = TDomain::dim;
1043 typedef typename TDomain::grid_type TGrid;
1044 typedef typename grid_dim_traits<dim>::element_type TElem;
1045 typedef typename grid_dim_traits<dim>::side_type TSide;
1046 typedef typename contrained_dim_traits<dim>::contrained_side_type TContrainedSide;
1047 typedef typename contrained_dim_traits<dim>::contraining_side_type TContrainingSide;
1048
1049 typedef typename TDomain::position_accessor_type TPositionAccessor;
1050
1051 const size_t fct = 0;
1052
1053 // get domain
1054 SmartPtr<TDomain> domain = spF->domain();
1055 SmartPtr<typename TDomain::grid_type> grid = domain->grid();
1056 typename TDomain::position_accessor_type& aaPos = domain->position_accessor();
1057
1058 // get dof distribution
1059 SmartPtr<DoFDistribution> dd = spF->dd();
1060
1061 // get iterators
1062 typename DoFDistribution::traits<TElem>::iterator iter, iterEnd;
1063 iter = dd->begin<TElem>(SurfaceView::ALL); // SurfaceView::MG_ALL
1064 iterEnd = dd->end<TElem>(SurfaceView::ALL);
1065
1066 // resize result vectors
1067 const size_t numIndex = spF->num_dofs();
1068 vStencil.resize(numIndex);
1069 vIndex.resize(numIndex);
1070
1071
1072 // storage (for reuse)
1073 std::vector<TElem*> vElem;
1074 std::vector<size_t> vNbrIndex;
1075 std::vector< MathVector<dim> > vDistance;
1076
1077
1079/*
1080 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spCond_F = spF->clone();
1081 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spCond_1 = spF->clone();
1082 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spCond_infty = spF->clone();
1083 DenseMatrix<VariableArray2<number> > Block;
1084*/
1086
1087
1088 DenseMatrix<VariableArray2<number> > BlockInv;
1089
1090 // loop all vertices
1091 for(;iter != iterEnd; ++iter)
1092 {
1093 // get vertex
1094 TElem* elem = *iter;
1095
1097 // Collect suitable neighbors for stencil creation
1099
1100 CollectStencilNeighbors_NeumannZeroBND_IndexAndDistance(vElem, vNbrIndex, vDistance, elem, *grid, aaPos, spF, fillElemSizeIntoVector);
1101
1102 const int numNeighborsRequired = 2*dim + (dim * (dim-1)) / 2;
1103 const int numSideNeighbors = std::pow(2,dim);
1104 if(vNbrIndex.size() < numNeighborsRequired || vDistance.size() < numNeighborsRequired)
1105 UG_THROW("Wrong number of neighbors detected: " << vNbrIndex.size());
1106
1107 // select additional neighbors beside the side connecting elems
1108 if(dim == 3) UG_THROW("DamageFunctionUpdater: This is 2d only --- extend to 3d by looking for 3 additional neighbors");
1109 if(numNeighborsRequired != numSideNeighbors+1) UG_THROW("DamageFunctionUpdater: Only 2d implemented")
1110
1111 // loop all "additional" neighbors (i.e. not those for sides on full refined meshes)
1112 // to find the clostest to choose
1113 int closest = -1; number closestDist = std::numeric_limits<number>::max();
1114 for(int k = numSideNeighbors; k < (int)vDistance.size(); ++k)
1115 {
1116 number dist = VecTwoNorm(vDistance[k]);
1117 if(dist < closestDist){closest = k; closestDist = dist;}
1118 }
1119 if(closest < 0) UG_THROW("DamageFunctionUpdater: closest not detected.")
1120
1121 vElem[numSideNeighbors] = vElem[closest];
1122 vNbrIndex[numSideNeighbors] = vNbrIndex[closest];
1123 vDistance[numSideNeighbors] = vDistance[closest];
1124
1126 // Create interpolation matrix and invert
1128
1129 BlockInv.resize(numNeighborsRequired, numNeighborsRequired);
1130 BlockInv = 0.0;
1131 for (size_t j = 0; j < numNeighborsRequired; ++j)
1132 {
1133 for (int d = 0; d < dim; ++d)
1134 BlockInv(j,d) = vDistance[j][d];
1135
1136 int cnt = dim;
1137 for (int d1 = 0; d1 < dim-1; ++d1)
1138 for (int d2 = d1+1; d2 < dim; ++d2)
1139 {
1140 BlockInv(j,cnt++) = vDistance[j][d1] * vDistance[j][d2];
1141 }
1142
1143 for (int d = 0; d < dim; ++d)
1144 BlockInv(j,cnt++) = 0.5 * vDistance[j][d] * vDistance[j][d];
1145
1146 if(cnt != numNeighborsRequired)
1147 UG_THROW("Wrong number of equations")
1148 }
1149
1151// Block = BlockInv;
1153
1154 if(!Invert(BlockInv))
1155 UG_THROW("Cannot invert block");
1156
1158/*
1159 for (size_t i = 0; i < numNeighborsRequired; ++i){
1160 for (size_t j = 0; j < numNeighborsRequired; ++j){
1161
1162 number res = 0.0;
1163 for(size_t k = 0; k < numNeighborsRequired; ++k)
1164 res += Block(i,k) * BlockInv(k,j);
1165
1166 if(i == j){
1167 if( fabs(res - 1.0) > 1e-10 )
1168 UG_THROW("Inverse wrong, should be 1.0, but is: " << res);
1169 } else {
1170 if( fabs(res - 0.0) > 1e-10 )
1171 UG_THROW("Inverse wrong, should be 0.0, but is: " << res);
1172 }
1173 }
1174 }
1175
1176 const number D_F = MatFrobeniusNorm(Block);
1177 const number DInv_F = MatFrobeniusNorm(BlockInv);
1178
1179 const number D_1 = MatOneNorm(Block);
1180 const number DInv_1 = MatOneNorm(BlockInv);
1181
1182 const number D_infty = MatInftyNorm(Block);
1183 const number DInv_infty = MatInftyNorm(BlockInv);
1184
1185 (*spCond_F)[i] = D_F * DInv_F;
1186 (*spCond_1)[i] = D_1 * DInv_1;
1187 (*spCond_infty)[i] = D_infty * DInv_infty;
1188*/
1190
1191
1192
1194 // extract second-order derivative subblock
1196
1197
1198 // add index
1199 std::vector<DoFIndex> ind;
1200 if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
1201 const size_t i = ind[0][0];
1202
1203 vStencil[i].clear();
1204 vStencil[i].resize(numNeighborsRequired+1);
1205 vStencil[i][0] = 0.0;
1206
1207 for (size_t k = 0; k < numNeighborsRequired; ++k){
1208 vStencil[i][k+1] = 0.0;
1209 for (int d = 0; d < dim; ++d){
1210 vStencil[i][k+1] += BlockInv( (numNeighborsRequired-dim)+d, k);
1211 vStencil[i][0] -= BlockInv( (numNeighborsRequired-dim)+d, k);
1212 }
1213 }
1214
1215 vIndex[i].clear();
1216 vIndex[i].resize(numNeighborsRequired+1);
1217 vIndex[i][0] = i;
1218 for (size_t k = 0; k < numNeighborsRequired; ++k){
1219 vIndex[i][k+1] = vNbrIndex[k];
1220 }
1221 }
1222
1223
1225/*
1226 write_debug(spCond_F, "Cond_D_F", 0, 0);
1227 write_debug(spCond_1, "Cond_D_1", 0, 0);
1228 write_debug(spCond_infty, "Cond_D_infty", 0, 0);
1229*/
1231
1232}
1233
1234
1235
1237// By least squares
1239
1240
1241template <typename TDomain>
1242void InitLaplacian_LeastSquares(
1243 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF,
1244 std::vector< std::vector< number > >& vStencil,
1245 std::vector< std::vector<size_t> >& vIndex, bool fillElemSizeIntoVector)
1246{
1247 PROFILE_BEGIN_GROUP(DamageFunctionUpdater_init_LeastSquares, "Small Strain Mech");
1248
1249 static const int dim = TDomain::dim;
1250 typedef typename TDomain::grid_type TGrid;
1251 typedef typename grid_dim_traits<dim>::element_type TElem;
1252 typedef typename grid_dim_traits<dim>::side_type TSide;
1253 typedef typename contrained_dim_traits<dim>::contrained_side_type TContrainedSide;
1254 typedef typename contrained_dim_traits<dim>::contraining_side_type TContrainingSide;
1255
1256 typedef typename TDomain::position_accessor_type TPositionAccessor;
1257
1258 const size_t fct = 0;
1259
1260 // get domain
1261 SmartPtr<TDomain> domain = spF->domain();
1262 SmartPtr<typename TDomain::grid_type> grid = domain->grid();
1263 typename TDomain::position_accessor_type& aaPos = domain->position_accessor();
1264
1265 // get dof distribution
1266 SmartPtr<DoFDistribution> dd = spF->dd();
1267
1268 // get iterators
1269 typename DoFDistribution::traits<TElem>::iterator iter, iterEnd;
1270 iter = dd->begin<TElem>(SurfaceView::ALL); // SurfaceView::MG_ALL
1271 iterEnd = dd->end<TElem>(SurfaceView::ALL);
1272
1273 // resize result vectors
1274 const size_t numIndex = spF->num_dofs();
1275 vStencil.resize(numIndex);
1276 vIndex.resize(numIndex);
1277
1278
1279 // storage (for reuse)
1280 std::vector<TElem*> vElem;
1281 std::vector<size_t> vNbrIndex;
1282 std::vector< MathVector<dim> > vDistance;
1283
1284 DenseMatrix<VariableArray2<number> > X, XTX, B;
1285 std::vector<number> W;
1286
1287 // loop all vertices
1288 for(;iter != iterEnd; ++iter)
1289 {
1290 // get vertex
1291 TElem* elem = *iter;
1292
1294 // Collect suitable neighbors for stencil creation
1296
1297 CollectStencilNeighbors_NeumannZeroBND_IndexAndDistance(vElem, vNbrIndex, vDistance, elem, *grid, aaPos, spF, fillElemSizeIntoVector);
1298
1299 const int numDerivs = 2*dim + (dim * (dim-1)) / 2; // or ( dim * (dim+3) ) / 2
1300 const int numNeighbors = vElem.size();
1301 if(vElem.size() < numDerivs || vNbrIndex.size() < numDerivs || vDistance.size() < numDerivs)
1302 UG_THROW("Wrong number of neighbors detected: #elems: " << vElem.size() <<
1303 ", #Index: " << vNbrIndex.size() << ", #Distance: " << vDistance.size());
1304
1306 // Create interpolation matrix and invert
1308
1309 // set weights
1310 W.resize(numNeighbors);
1311 for (int i = 0; i < numNeighbors; ++i){
1312// W[i] = 1.0;
1313// W[i] = 1.0/VecTwoNorm(vDistance[i]);
1314 W[i] = 1.0/VecTwoNormSq(vDistance[i]);
1315 }
1316
1317 // set X
1318 X.resize(numNeighbors, numDerivs);
1319 X = 0.0;
1320 for (int j = 0; j < numNeighbors; ++j)
1321 {
1322 for (int d = 0; d < dim; ++d)
1323 X(j,d) = vDistance[j][d];
1324
1325 int cnt = dim;
1326 for (int d1 = 0; d1 < dim-1; ++d1)
1327 for (int d2 = d1+1; d2 < dim; ++d2)
1328 {
1329 X(j,cnt++) = vDistance[j][d1] * vDistance[j][d2];
1330 }
1331
1332 for (int d = 0; d < dim; ++d)
1333 X(j,cnt++) = 0.5 * vDistance[j][d] * vDistance[j][d];
1334
1335 if(cnt != numDerivs)
1336 UG_THROW("Wrong number of equations")
1337 }
1338
1339 // W * X
1340 for (int i = 0; i < numNeighbors; ++i){
1341 for(int j = 0; j < numDerivs; ++j){
1342 X(i,j) *= W[i];
1343 }
1344 }
1345
1346 // XTX = (WX)^T * WX = (X^T W^T W X)
1347 XTX.resize(numDerivs, numDerivs);
1348 for(int i = 0; i < numDerivs; ++i){
1349 for(int j = 0; j < i; ++j){
1350 XTX(i,j) = 0;
1351 for(int k = 0; k < numNeighbors; ++k){
1352 XTX(i,j) += X(k,i) * X(k,j);
1353 }
1354 XTX(j,i) = XTX(i,j);
1355 }
1356 XTX(i,i) = 0;
1357 for(int k = 0; k < numNeighbors; ++k){
1358 XTX(i,i) += X(k,i) * X(k,i);
1359 }
1360 }
1361
1362 // XTX = (X^T W^T W X)^{-1}
1363 if(!Invert(XTX))
1364 UG_THROW("Cannot invert block");
1365
1366 // B = (X^T W^T W X)^{-1} * X^T * W^T * W
1367 B.resize(numDerivs, numNeighbors);
1368 for(int i = 0; i < numDerivs; ++i){
1369 for(int j = 0; j < numNeighbors; ++j){
1370 B(i,j) = 0;
1371 for(int k = 0; k < numDerivs; ++k){
1372 B(i,j) += XTX(i,k) * X(j,k);
1373 }
1374 B(i,j) *= W[j];
1375 }
1376 }
1377
1378
1380 // extract second-order derivative subblock
1382
1383 // add index
1384 std::vector<DoFIndex> ind;
1385 if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
1386 const size_t i = ind[0][0];
1387
1388 vStencil[i].clear();
1389 vStencil[i].resize(numNeighbors+1);
1390 vStencil[i][0] = 0.0;
1391
1392 for (int k = 0; k < numNeighbors; ++k){
1393 vStencil[i][k+1] = 0.0;
1394 for (int d = 0; d < dim; ++d){
1395 vStencil[i][k+1] += B( (numDerivs-dim)+d, k);
1396 vStencil[i][0] -= B( (numDerivs-dim)+d, k);
1397 }
1398 }
1399
1400 vIndex[i].clear();
1401 vIndex[i].resize(numNeighbors+1);
1402 vIndex[i][0] = i;
1403 for (int k = 0; k < numNeighbors; ++k){
1404 vIndex[i][k+1] = vNbrIndex[k];
1405 }
1406 }
1407
1408}
1409
1410
1412// By Taylor for discrete laplacian
1414
1415
1416template <typename TDomain>
1417void InitLaplacian_TaylorDirect(
1418 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF,
1419 std::vector< std::vector< number > >& vStencil,
1420 std::vector< std::vector<size_t> >& vIndex, bool fillElemSizeIntoVector = false)
1421{
1422 PROFILE_BEGIN_GROUP(DamageFunctionUpdater_init_TaylorDirect, "Small Strain Mech");
1423
1424 static const int dim = TDomain::dim;
1425 typedef typename TDomain::grid_type TGrid;
1426 typedef typename grid_dim_traits<dim>::element_type TElem;
1427 typedef typename grid_dim_traits<dim>::side_type TSide;
1428 typedef typename contrained_dim_traits<dim>::contrained_side_type TContrainedSide;
1429 typedef typename contrained_dim_traits<dim>::contraining_side_type TContrainingSide;
1430
1431 typedef typename TDomain::position_accessor_type TPositionAccessor;
1432
1433 const size_t fct = 0;
1434
1435 // get domain
1436 SmartPtr<TDomain> domain = spF->domain();
1437 SmartPtr<typename TDomain::grid_type> grid = domain->grid();
1438 typename TDomain::position_accessor_type& aaPos = domain->position_accessor();
1439
1440 // get dof distribution
1441 SmartPtr<DoFDistribution> dd = spF->dd();
1442
1443 // get iterators
1444 typename DoFDistribution::traits<TElem>::iterator iter, iterEnd;
1445 iter = dd->begin<TElem>(SurfaceView::ALL); // SurfaceView::MG_ALL
1446 iterEnd = dd->end<TElem>(SurfaceView::ALL);
1447
1448 // resize result vectors
1449 const size_t numIndex = spF->num_dofs();
1450 vStencil.resize(numIndex);
1451 vIndex.resize(numIndex);
1452
1453
1454 // storage (for reuse)
1455 std::vector<TElem*> vElem;
1456 std::vector<size_t> vNbrIndex;
1457 std::vector< MathVector<dim> > vDistance;
1458
1459 DenseMatrix<VariableArray2<number> > X, XTX;
1460 std::vector<number> XTXr;
1461
1462 // loop all vertices
1463 for(;iter != iterEnd; ++iter)
1464 {
1465 // get vertex
1466 TElem* elem = *iter;
1467
1469 // Collect suitable neighbors for stencil creation
1471
1472 CollectStencilNeighbors_NeumannZeroBND_IndexAndDistance(vElem, vNbrIndex, vDistance, elem, *grid, aaPos, spF, fillElemSizeIntoVector);
1473
1474 const int numDerivs = 2*dim + (dim * (dim-1)) / 2; // or ( dim * (dim+3) ) / 2
1475 const int numNeighbors = vElem.size();
1476 if(vElem.size() < numDerivs || vNbrIndex.size() < numDerivs || vDistance.size() < numDerivs)
1477 UG_THROW("Wrong number of neighbors detected: #elems: " << vElem.size() <<
1478 ", #Index: " << vNbrIndex.size() << ", #Distance: " << vDistance.size());
1479
1481 // Create interpolation matrix and invert
1483
1484 // set X
1485 X.resize(numNeighbors, numDerivs);
1486 X = 0.0;
1487 for (int j = 0; j < numNeighbors; ++j)
1488 {
1489 for (int d = 0; d < dim; ++d)
1490 X(j,d) = vDistance[j][d];
1491
1492 int cnt = dim;
1493 for (int d1 = 0; d1 < dim-1; ++d1)
1494 for (int d2 = d1+1; d2 < dim; ++d2)
1495 {
1496 X(j,cnt++) = vDistance[j][d1] * vDistance[j][d2];
1497 }
1498
1499 for (int d = 0; d < dim; ++d)
1500 X(j,cnt++) = 0.5 * vDistance[j][d] * vDistance[j][d];
1501
1502 if(cnt != numDerivs)
1503 UG_THROW("Wrong number of equations")
1504 }
1505
1506 // X^T X
1507 XTX.resize(numDerivs, numDerivs);
1508 for(int i = 0; i < numDerivs; ++i){
1509 for(int j = 0; j < i; ++j){
1510 XTX(i,j) = 0;
1511 for(int k = 0; k < numNeighbors; ++k){
1512 XTX(i,j) += X(k,i) * X(k,j);
1513 }
1514 XTX(j,i) = XTX(i,j);
1515 }
1516 XTX(i,i) = 0;
1517 for(int k = 0; k < numNeighbors; ++k){
1518 XTX(i,i) += X(k,i) * X(k,i);
1519 }
1520 }
1521
1522 // (X^T X)^{-1}
1523 if(!Invert(XTX))
1524 UG_THROW("Cannot invert block");
1525
1527 // extract second-order derivative subblock
1529
1530 // add index
1531 std::vector<DoFIndex> ind;
1532 if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
1533 const size_t i = ind[0][0];
1534
1535 vStencil[i].clear();
1536 vStencil[i].resize(numNeighbors+1);
1537 vStencil[i][0] = 0.0;
1538
1539 // XTXr = (X^T X)^{-1} * s
1540 XTXr.resize(numNeighbors);
1541 for (int j = 0; j < numNeighbors; ++j){
1542
1543 XTXr[j] = 0.0;
1544 for (int d = 0; d < dim; ++d)
1545 XTXr[j] += XTX(j,(numDerivs-dim)+d);
1546 }
1547
1548 // X * (X^T X)^{-1} * s
1549 for (int k = 0; k < numNeighbors; ++k){
1550 vStencil[i][k+1] = 0.0;
1551 for (int j = 0; j < numNeighbors; ++j){
1552 vStencil[i][k+1] += X(k,j) * XTXr[j];
1553 }
1554 vStencil[i][0] -= vStencil[i][k+1];
1555 }
1556
1557 vIndex[i].clear();
1558 vIndex[i].resize(numNeighbors+1);
1559 vIndex[i][0] = i;
1560 for (int k = 0; k < numNeighbors; ++k){
1561 vIndex[i][k+1] = vNbrIndex[k];
1562 }
1563 }
1564
1565}
1566
1567
1568
1572// DamageFunctionUpdater
1576
1577
1578template <typename TDomain>
1579void DamageFunctionUpdater<TDomain>::
1580set_disc_type(const std::string& theType)
1581{
1582 std::string type = TrimString(theType);
1583 std::transform(type.begin(), type.end(), type.begin(), ::tolower);
1584
1585 UG_LOG ("DamageFunctionUpdater: setting laplacian type to '"<< theType <<"' \n");
1586 if(type == "least-squares") {m_discType = _LEAST_SQUARES_; return;}
1587 if(type == "taylor-expansion") {m_discType = _TAYLOR_EXPANSION_; return;}
1588 if(type == "partial-integration") {m_discType = _PARTIAL_INTEGRATION_; return;}
1589 if(type == "taylor-direct") {m_discType = _TAYLOR_DIRECT_; return;}
1590
1591 UG_THROW("DamageFunctionUpdater: unrecognized type '"<<type<<"'.");
1592}
1593
1594
1595template <typename TDomain>
1596bool DamageFunctionUpdater<TDomain>::
1597solve( SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF,
1598 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spPsi0,
1599 const number beta, const number r,
1600 const number eps, const int maxIter, const number dampNewton)
1601{
1602 PROFILE_BEGIN_GROUP(DamageFunctionUpdater_solve, "Small Strain Mech");
1603
1605 // check if has to be rebuild
1607
1608 static int call = 0; call++;
1609
1610 // get approximation space
1611 ConstSmartPtr<ApproximationSpace<TDomain> > approxSpace = spF->approx_space();
1612 if(approxSpace != spPsi0->approx_space())
1613 UG_THROW("DamageFunctionUpdater<TDomain>::solve: expected same ApproximationSpace for f and psi0");
1614
1615 // check revision counter if grid / approx space has changed since last call
1616 if(m_ApproxSpaceRevision != approxSpace->revision())
1617 {
1618 UG_LOG ("DamageFunctionUpdater: reinit laplacian using type '"<< m_discType <<"'\n");
1619 PROFILE_BEGIN_GROUP(DamageFunctionUpdater_init, "Small Strain Mech");
1620 // (re-)initialize setting
1621 switch(m_discType){
1622 case _PARTIAL_INTEGRATION_:
1623 InitLaplacian_PartialIntegration(spF, m_vStencil, m_vIndex, m_quadRuleType);
1624 break;
1625 case _TAYLOR_EXPANSION_:
1626 InitLaplacian_TaylorExpansion(spF, m_vStencil, m_vIndex);
1627 break;
1628 case _LEAST_SQUARES_:
1629 InitLaplacian_LeastSquares(spF, m_vStencil, m_vIndex);
1630 break;
1631 case _TAYLOR_DIRECT_:
1632 InitLaplacian_TaylorDirect(spF, m_vStencil, m_vIndex);
1633 break;
1634 default:
1635 UG_THROW("DamageFunctionUpdater: internal error, "
1636 "unrecognized type number '"<<m_discType<<"'.");
1637 }
1638
1639 // remember revision counter of approx space
1640 m_ApproxSpaceRevision = approxSpace->revision();
1641
1642 write_stencil_matrix_debug(spF, "STENCIL", call);
1643 }
1644
1645
1646
1648 // apply newton method
1650
1651// const size_t numElem = m_vIndex.size();
1652// const number sqrtNumElem = sqrt(numElem);
1653
1654 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spLambdaOld = spF->clone();
1655 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spPhi = spF->clone();
1656
1657
1658 // number normPhi = std::numeric_limits<number>::max();
1659 int iterCnt = 0;
1660
1662/* for(size_t i = 0; i < m_vIndex.size(); ++i)
1663 (*spLambdaOld)[i] = DLambda(i);
1664 write_debug(spLambdaOld, "DLambda", call, iterCnt);
1665
1666 write_debug(spPsi0, "Psi0", call, iterCnt);
1667 write_debug(spF, "F", call, iterCnt);
1668*/
1669
1670 number maxPhi = std::numeric_limits<number>::max();
1671
1672// while(normPhi > eps * sqrtNumElem && (iterCnt++ <= maxIter) )
1673 while(maxPhi > eps && (iterCnt++ <= maxIter) )
1674 {
1675 // normPhi = 0.0;
1676 maxPhi = 0.0;
1677
1678 for(size_t i = 0; i < m_vIndex.size(); ++i)
1679 (*spLambdaOld)[i] = Lambda(i, spF);
1680
1682// write_debug(spLambdaOld, "Lambda", call, iterCnt);
1684
1685 for(size_t i = 0; i < m_vIndex.size(); ++i)
1686 {
1687 const number lambda = (*spLambdaOld)[i];
1688 const number& psi0 = (*spPsi0)[i];
1689 number& f = (*spF)[i];
1690
1691 number phi = f * ( psi0 - beta * lambda) - r;
1692
1694// (*spPhi)[i] = phi;
1696
1697
1698 if(phi < eps){
1699// phi = 0;
1700// normPhi += 0.0 * 0.0;
1701 } else {
1702
1703 // ORIGINAL
1704 // f = f - (phi / (psi0 - beta * (lambda + f * DLambda(i)) ));
1705
1706
1707 // DAMPED NEWTON
1708 f = f - dampNewton * (phi / (psi0 - beta * (lambda + f * DLambda(i)) ));
1709
1710 // QUASI-NEWTON
1711 //f = f - (1/10)*(phi / (psi0 - beta * lambda));
1712
1713 // FIXPOINT
1714 //f = phiScale * (f * ( psi0 - beta * lambda) - r) + f;
1715
1716
1718 if(std::isfinite(f) == false){
1719
1720 UG_LOG(" ############### \n");
1721 UG_LOG("f : " << f << "\n");
1722 UG_LOG("psi0 : " << psi0 << "\n");
1723 UG_LOG("lambda: " << lambda << "\n");
1724 UG_LOG("DLambda: " << DLambda(i) << "\n");
1725 UG_LOG("beta : " << beta << "\n");
1726 UG_LOG("r : " << r << "\n");
1727 UG_LOG("phi : " << phi << "\n");
1728 UG_LOG(" ############### \n");
1729 UG_THROW("Value for f not finite, but: " << f);
1730 }
1732
1733 // normPhi += phi*phi;
1734 maxPhi = std::max(maxPhi, phi);
1735 }
1736 }
1737
1739// write_debug(spPhi, "Phi", call, iterCnt);
1740// write_debug(spF, "F", call, iterCnt);
1742
1743 // normPhi = sqrt(normPhi);
1744// UG_LOG(" ######################### (end sweep) ################################### \n");
1745// UG_LOG ("DamageFunctionUpdater: normPhi: " << normPhi << "\n");
1746 }
1747
1748 m_lastNumIters = iterCnt;
1749
1750 if(iterCnt >= maxIter){
1751 // UG_THROW("DamageFunctionUpdater: no convergence after " << iterCnt << " iterations");
1752 UG_LOG("DamageFunctionUpdater: no convergence after " << iterCnt << " iterations");
1753 return false;
1754 }
1755
1756
1757 UG_LOG ("DamageFunctionUpdater: maxPhi: " << maxPhi << " after " <<iterCnt << " iterations\n");
1758
1759 return true;
1760}
1761
1762
1763
1764
1765template <typename TDomain>
1766void DamageFunctionUpdater<TDomain>::
1767set_debug(SmartPtr<GridFunctionDebugWriter<TDomain, CPUAlgebra> > spDebugWriter)
1768{
1769 m_spDebugWriter = spDebugWriter;
1770}
1771
1772template <typename TDomain>
1773void DamageFunctionUpdater<TDomain>::
1774write_debug(SmartPtr<GridFunction<TDomain, CPUAlgebra> > spGF, std::string name,
1775 int call, int iter)
1776{
1777 if(m_spDebugWriter.invalid()) return;
1778
1779// build name
1780 GridLevel gl = spGF->grid_level();
1781 std::stringstream ss;
1782 ss << "InDamageUpdate_" << name ;
1783 ss << "_call" << std::setfill('0') << std::setw(3) << call;
1784 if (iter >= 0) ss << "_iter" << std::setfill('0') << std::setw(3) << iter;
1785 ss << ".vec";
1786
1787// write
1788 m_spDebugWriter->set_grid_level(gl);
1789 m_spDebugWriter->write_vector(*spGF, ss.str().c_str());
1790}
1791
1792template <typename TDomain>
1793void DamageFunctionUpdater<TDomain>::
1794write_stencil_matrix_debug(
1795 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spGF,
1796 std::string name, int call)
1797{
1798 if(m_spDebugWriter.invalid()) return;
1799
1800// build name
1801 GridLevel gl = spGF->grid_level();
1802 int numDoFs = spGF->num_dofs();
1803 std::stringstream ss;
1804 ss << "InDamageUpdate_" << name ;
1805 ss << "_call" << std::setfill('0') << std::setw(3) << call;
1806 //if (iter >= 0) ss << "_iter" << std::setfill('0') << std::setw(3) << iter;
1807 ss << ".mat";
1808
1809 typedef CPUAlgebra::matrix_type TMat;
1810
1811 TMat A;
1812 A.resize_and_clear(numDoFs,numDoFs);
1813
1814 for(size_t i = 0; i < m_vStencil.size(); ++i){
1815 for(size_t k = 0; k < m_vStencil[i].size(); ++k){
1816 const int j = m_vIndex[i][k];
1817
1818 A(i,j) = m_vStencil[i][k];
1819 }
1820 }
1821
1822// write
1823 m_spDebugWriter->set_grid_level(gl);
1824 m_spDebugWriter->write_matrix(A, ss.str().c_str());
1825
1826}
1827
1828
1829
1833// RelativeDensityUpdater
1837
1838template <typename TDomain>
1839void RelativeDensityUpdater<TDomain>::
1840set_disc_type(const std::string& theType)
1841{
1842 std::string type = TrimString(theType);
1843 std::transform(type.begin(), type.end(), type.begin(), ::tolower);
1844
1845 UG_LOG ("DamageFunctionUpdater: setting laplacian type to '"<< theType <<"' \n");
1846 if(type == "least-squares") {m_discType = _LEAST_SQUARES_; return;}
1847 if(type == "taylor-expansion") {m_discType = _TAYLOR_EXPANSION_; return;}
1848 if(type == "partial-integration") {m_discType = _PARTIAL_INTEGRATION_; return;}
1849 if(type == "taylor-direct") {m_discType = _TAYLOR_DIRECT_; return;}
1850
1851 UG_THROW("DamageFunctionUpdater: unrecognized type '"<<type<<"'.");
1852}
1853
1854
1855
1856template <typename TDomain>
1857std::vector<number> RelativeDensityUpdater<TDomain>::
1858solve( SmartPtr<GridFunction<TDomain, CPUAlgebra> > spChi,
1859 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spDrivingForce,
1860 const number betaStar, const number etaChiStar,
1861 const number chiMin, const number dt, const int p,
1862 const number rho_target, const number MassTol)
1863{
1864 PROFILE_BEGIN_GROUP(RelativeDensityUpdater_solve, "Small Strain Mech");
1865 static int call = 0; call++;
1866
1867
1868
1870 // rebuild after grid adaption
1872
1873 // get approximation space
1874 ConstSmartPtr<ApproximationSpace<TDomain> > approxSpace = spChi->approx_space();
1875 if(approxSpace != spDrivingForce->approx_space())
1876 UG_THROW("RelativeDensityUpdater<TDomain>::solve: expected same ApproximationSpace for f and psi0");
1877
1878 // check revision counter if grid / approx space has changed since last call
1879 if(m_ApproxSpaceRevision != approxSpace->revision())
1880 {
1881 m_spElemSize = spChi->clone();
1882 m_spLaplaceChi = spChi->clone();
1883 m_spChiTrial = spChi->clone();
1884
1885
1886 UG_LOG ("RelativeDensityUpdater: reinit laplacian using type '"<< m_discType <<"'\n");
1887 PROFILE_BEGIN_GROUP(DamageFunctionUpdater_init, "Small Strain Mech");
1888 // (re-)initialize setting
1889 switch(m_discType){
1890 case _PARTIAL_INTEGRATION_:
1891 InitLaplacian_PartialIntegration(m_spElemSize, m_vStencil, m_vIndex, m_quadRuleType, true);
1892 break;
1893 case _TAYLOR_EXPANSION_:
1894 InitLaplacian_TaylorExpansion(m_spElemSize, m_vStencil, m_vIndex, true);
1895 break;
1896 case _LEAST_SQUARES_:
1897 InitLaplacian_LeastSquares(m_spElemSize, m_vStencil, m_vIndex, true);
1898 break;
1899 case _TAYLOR_DIRECT_:
1900 InitLaplacian_TaylorDirect(m_spElemSize, m_vStencil, m_vIndex, true);
1901 break;
1902 default:
1903 UG_THROW("RelativeDensityUpdater: internal error, "
1904 "unrecognized type number '"<<m_discType<<"'.");
1905 }
1906
1907 // remember revision counter of approx space
1908 m_ApproxSpaceRevision = approxSpace->revision();
1909
1910 // debug output
1911 write_stencil_matrix_debug(spChi, "STENCIL", call);
1912 }
1913
1914
1916 // loop until maxiter
1918
1919 // find smallest h^2
1920 number h2 = std::numeric_limits<number>::max();
1921 for(size_t i = 0; i < m_vIndex.size(); ++i){
1922 const number vol = (*m_spElemSize)[i];
1923 h2 = std::min(h2, sqrt(dim) * std::pow(vol, 2.0/dim));
1924 }
1925
1926 const int n = std::floor(6*betaStar / (etaChiStar * h2)) + 1;
1927
1928 const number dt_chi = dt / n;
1929
1930 int numBisect = 0;
1931
1932 // loop over inner (smaller) timesteps
1933 for(int j = 0; j < n; ++j)
1934 {
1936 // compute parameter beta, eta
1938
1939 number min_p_chi = std::numeric_limits<number>::max();
1940 number max_p_chi = std::numeric_limits<number>::min();
1941
1942 number int_g_p = 0.0, int_g = 0.0;
1943 for(size_t i = 0; i < m_vIndex.size(); ++i)
1944 {
1945 const number chi = (*spChi)[i];
1946 const number vol = (*m_spElemSize)[i];
1947 const number p_chi = (*spDrivingForce)[i];
1948 (*m_spLaplaceChi)[i] = Lambda(i, spChi);
1949
1950 const number g = (chi - chiMin) * (1.0 - chi);
1951
1952 int_g += vol * g;
1953 int_g_p += vol * g * p_chi;
1954
1955 min_p_chi = std::min(min_p_chi, p_chi);
1956 max_p_chi = std::max(max_p_chi, p_chi);
1957 }
1958
1959 const number p_w = int_g_p / int_g;
1960
1961 number beta = betaStar * p_w;
1962 const number eta_chi = etaChiStar * p_w;
1963
1964
1966 // compute lagrange multiplier lambda
1968/*
1969 number int_1 = 0.0, int_partPsi = 0.0, int_LaplaceChi = 0.0;
1970 for(size_t i = 0; i < m_vIndex.size(); ++i)
1971 {
1972 const number vol = (*m_spElemSize)[i];
1973 const number p_chi = (*spDrivingForce)[i];
1974 const number laplaceChi = (*m_spLaplaceChi)[i] = Lambda(i, spChi);
1975
1976 int_1 += vol;
1977 int_partPsi += p_chi * vol;
1978 int_LaplaceChi += laplaceChi * vol;
1979 }
1980
1981 const number lambda = (1.0 / int_1) * (int_partPsi + beta * int_LaplaceChi);
1982*/
1984 // update relative density (and theredy also driving force)
1986
1987
1989// write_debug(m_spLaplaceChi, "LaplaceChi", call, -1);
1991 number Llow = (min_p_chi - eta_chi); // untere Grenze für Lagrangemultiplikator
1992 number Lhigh = (max_p_chi + eta_chi); // obere Grenze für Lagrangemultiplikator
1993 number L_tr = 0.0; // why not: (Lhigh + Llow)/2.0; ???
1994
1995 // bisection to find lagrange param
1996 while(true)
1997 {
1998 numBisect++;
1999
2000 number Vol_Omega = 0.0;
2001 number rho_tr = 0.0;
2002
2003 // Evolution equation for each element (trial step)
2004 for(size_t i = 0; i < m_vIndex.size(); ++i)
2005 {
2006 const number chi = (*spChi)[i];
2007 const number p_chi = (*spDrivingForce)[i];
2008 const number laplaceChi = (*m_spLaplaceChi)[i];
2009 const number vol = (*m_spElemSize)[i];
2010
2011 if(m_bEnforceLocalRequiredBeta){
2012 beta = std::max(2 * sqrt(dim) * std::pow(vol, 2.0/dim) * p_w, betaStar * p_w);
2013 }
2014
2015 number& chi_tr = (*m_spChiTrial)[i];
2016
2017 // update relative density
2018 chi_tr = chi + (dt_chi / eta_chi) * (p_chi - L_tr + beta * laplaceChi);
2019
2021 if(std::isfinite(chi_tr) == false){
2022
2023 UG_LOG(" ############### \n");
2024 UG_LOG("chi_tr : " << chi_tr << "\n");
2025 UG_LOG("dt_chi : " << dt_chi << "\n");
2026 UG_LOG("eta_chi : " << eta_chi << "\n");
2027 UG_LOG("p_chi : " << p_chi << "\n");
2028 UG_LOG("lambda : " << L_tr << "\n");
2029 UG_LOG("beta : " << beta << "\n");
2030 UG_LOG("laplaceChi : " << laplaceChi << "\n");
2031 UG_LOG("p_w : " << p_w << "\n");
2032 UG_LOG(" ############### \n");
2033 UG_THROW("Value for chi not finite, but: " << chi_tr);
2034 }
2036
2037 if(chi_tr > 1.0) chi_tr = 1.0;
2038 if(chi_tr < chiMin) chi_tr = chiMin;
2039
2040 rho_tr += chi_tr * vol;
2041 Vol_Omega += vol;
2042 }
2043
2044 rho_tr /= Vol_Omega;
2045
2046 // Volume constraint
2047 // trial Volume
2048 if(fabs(rho_tr - rho_target) < MassTol) break; // V_tr = V_target
2049
2050 // Lagrange update
2051 if(rho_tr > rho_target) Llow = L_tr;
2052 if(rho_tr < rho_target) Lhigh = L_tr;
2053
2054 L_tr = (Lhigh + Llow)/2.0;
2055 }
2056
2057
2058 // accept update
2059 for(size_t i = 0; i < m_vIndex.size(); ++i)
2060 {
2061 number& chi = (*spChi)[i];
2062 number& p_chi = (*spDrivingForce)[i];
2063 const number chi_tr = (*m_spChiTrial)[i];
2064
2065
2066 p_chi *= std::pow( (chi_tr / chi), p-1);
2067
2068 chi = chi_tr;
2069 }
2070 }
2071
2072
2073 std::vector<number> vRes;
2074
2075 vRes.push_back((number)n);
2076 vRes.push_back((number)numBisect);
2077 vRes.push_back((number)numBisect/(number)n);
2078
2079 return vRes;
2080}
2081
2082
2083
2084
2085template <typename TDomain>
2086void RelativeDensityUpdater<TDomain>::
2087set_debug(SmartPtr<GridFunctionDebugWriter<TDomain, CPUAlgebra> > spDebugWriter)
2088{
2089 m_spDebugWriter = spDebugWriter;
2090}
2091
2092template <typename TDomain>
2093void RelativeDensityUpdater<TDomain>::
2094write_debug(SmartPtr<GridFunction<TDomain, CPUAlgebra> > spGF, std::string name,
2095 int call, int iter)
2096{
2097 if(m_spDebugWriter.invalid()) return;
2098
2099// build name
2100 GridLevel gl = spGF->grid_level();
2101 std::stringstream ss;
2102 ss << "InDensityUpdate_" << name ;
2103 ss << "_call" << std::setfill('0') << std::setw(3) << call;
2104 if (iter >= 0) ss << "_iter" << std::setfill('0') << std::setw(3) << iter;
2105 ss << ".vec";
2106
2107// write
2108 m_spDebugWriter->set_grid_level(gl);
2109 m_spDebugWriter->write_vector(*spGF, ss.str().c_str());
2110}
2111
2112template <typename TDomain>
2113void RelativeDensityUpdater<TDomain>::
2114write_stencil_matrix_debug(
2115 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spGF,
2116 std::string name, int call)
2117{
2118 if(m_spDebugWriter.invalid()) return;
2119
2120// build name
2121 GridLevel gl = spGF->grid_level();
2122 int numDoFs = spGF->num_dofs();
2123 std::stringstream ss;
2124 ss << "InDensityUpdate_" << name ;
2125 ss << "_call" << std::setfill('0') << std::setw(3) << call;
2126 //if (iter >= 0) ss << "_iter" << std::setfill('0') << std::setw(3) << iter;
2127 ss << ".mat";
2128
2129 typedef CPUAlgebra::matrix_type TMat;
2130
2131 TMat A;
2132 A.resize_and_clear(numDoFs,numDoFs);
2133
2134 for(size_t i = 0; i < m_vStencil.size(); ++i){
2135 for(size_t k = 0; k < m_vStencil[i].size(); ++k){
2136 const int j = m_vIndex[i][k];
2137
2138 A(i,j) = m_vStencil[i][k];
2139 }
2140 }
2141
2142// write
2143 m_spDebugWriter->set_grid_level(gl);
2144 m_spDebugWriter->write_matrix(A, ss.str().c_str());
2145
2146}
2147
2148
2149
2150
2154// Marking
2158
2159template<typename TDomain>
2160void MarkForAdaption_ValueRangeIndicator(
2161 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spChi,
2162 IRefiner& refiner,
2163 number lowerValueToCoarsen,
2164 number minValueToRefine, number maxValueToRefine,
2165 number upperValueToCoarsen,
2166 number maxJumpDifferenceToNeighborsForCoarsen,
2167 number minJumpDifferenceToNeighborsForRefine,
2168 int maxLevel)
2169{
2170 PROFILE_FUNC_GROUP("Small Strain Mech");
2171
2172 static const int dim = TDomain::dim;
2173 typedef typename grid_dim_traits<dim>::element_type TElem;
2174 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
2175 typedef typename TDomain::position_accessor_type position_accessor_type;
2176 //position_accessor_type& aaPos = spChi->domain()->position_accessor();
2177
2178 const int fct = 0; // \todo: generalize
2179
2181 // checks
2183
2184 if( minValueToRefine > maxValueToRefine ) UG_THROW("minValueToRefine > maxValueToRefine")
2185 if( lowerValueToCoarsen > minValueToRefine ) UG_THROW("lowerValueToCoarsen > minValueToRefine")
2186 if( upperValueToCoarsen < maxValueToRefine ) UG_THROW("upperValueToCoarsen < maxValueToRefine")
2187
2188
2189
2190
2191 // marking
2193
2194 // reset counter
2195 int numMarkedRefine = 0, numMarkedCoarse = 0;
2196
2197
2198 std::vector<MathVector<dim> > vCornerCoords;
2199 MathVector<dim> ElemCenter;
2200
2201 std::vector<TElem*> vNeighbors;
2202
2203 // loop elements for marking
2204 const_iterator iter = spChi->template begin<TElem>();
2205 const_iterator iterEnd = spChi->template end<TElem>();
2206 for(; iter != iterEnd; ++iter)
2207 {
2208 // get element
2209 TElem* elem = *iter;
2210
2211 std::vector<DoFIndex> ind;
2212 if(spChi->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
2213 const number val = DoFRef(*spChi, ind[0]);
2214
2215
2216 CollectSurfaceNeighbors(spChi, elem, vNeighbors);
2217
2218 // check whether all children of e of type TElem have similar values
2219 bool neighborsHaveLargeJump = false;
2220 bool neighborsHaveSmallJump = true;
2221 for(size_t i = 0; i < vNeighbors.size(); ++i){
2222
2223 TElem* neighbor = vNeighbors[i];
2224
2225 std::vector<DoFIndex> ind;
2226 if(spChi->inner_dof_indices(neighbor, fct, ind) != 1) UG_THROW("Wrong number dofs");
2227 const number valNeighbor = DoFRef(*spChi, ind[0]);
2228
2229 if(fabs(val - valNeighbor) > minJumpDifferenceToNeighborsForRefine)
2230 {
2231 neighborsHaveLargeJump = true;
2232 }
2233
2234 if(fabs(val - valNeighbor) > maxJumpDifferenceToNeighborsForCoarsen)
2235 {
2236 neighborsHaveSmallJump = false;
2237 }
2238
2239 }
2240
2241 // marks for refinement
2242 if( (val >= minValueToRefine && val <= maxValueToRefine)
2243 || (neighborsHaveLargeJump) )
2244 {
2245
2246 if(spChi->dd()->multi_grid()->get_level(elem) < maxLevel)
2247 {
2248 refiner.mark(elem, RM_REFINE);
2249 numMarkedRefine++;
2250 }
2251
2252 continue;
2253 }
2254
2255 // marks for coarsening
2256 if( (val < lowerValueToCoarsen || val > upperValueToCoarsen)
2257 && neighborsHaveSmallJump )
2258 {
2259
2260 // get the parent
2261 TElem* parent = dynamic_cast<TElem*>(spChi->dd()->multi_grid()->get_parent(elem));
2262 if(parent){
2263
2264 refiner.mark(elem, RM_COARSEN);
2265 numMarkedCoarse++;
2266
2267 }
2268
2269 }
2270
2271 }
2272
2274 // print infos
2276 UG_LOG("MarkForAdaption_ValueRangeIndicator:\n")
2277#ifdef UG_PARALLEL
2278 if(pcl::NumProcs() > 1)
2279 {
2280 UG_LOG(" >>> Marked on Proc "<<pcl::ProcRank()<<": refine: " << numMarkedRefine << ", coarsen : " << numMarkedCoarse << "\n");
2281 pcl::ProcessCommunicator com;
2282 int numMarkedRefineLocal = numMarkedRefine, numMarkedCoarseLocal = numMarkedCoarse;
2283 numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
2284 numMarkedCoarse = com.allreduce(numMarkedCoarseLocal, PCL_RO_SUM);
2285 }
2286#endif
2287
2288 UG_LOG(" >>> Marked refine: " << numMarkedRefine << ", coarsen: " << numMarkedCoarse << "\n" );
2289}
2290
2291
2292template<typename TDomain>
2293void MarkDamage( SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF,
2294 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spPsi0,
2295 IRefiner& refiner,
2296 number minValueToRefine, number maxValueToCoarsen,
2297 int maxLevel,
2298 const std::vector<MathVector<TDomain::dim,number>* >& vCenter,
2299 const std::vector<number>& vRadius)
2300{
2301 PROFILE_FUNC_GROUP("Small Strain Mech");
2302
2303 static const int dim = TDomain::dim;
2304 typedef typename grid_dim_traits<dim>::element_type TElem;
2305 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
2306 typedef typename TDomain::position_accessor_type position_accessor_type;
2307 position_accessor_type& aaPos = spF->domain()->position_accessor();
2308
2309 const int fct = 0; // \todo: generalize
2310
2312 // checks
2314
2315 if( (minValueToRefine < maxValueToCoarsen) )
2316 UG_THROW("Coarsen threshold must be smaller than refine threshold")
2317
2318 if( !(vRadius.size() == vCenter.size()) )
2319 UG_THROW("Num radius and center must match")
2320
2321
2322
2323 // marking
2325
2326 // we'll compare against the square radius.
2327 std::vector<number> vRadiusSq;
2328 for(size_t i = 0; i < vRadius.size(); ++i)
2329 vRadiusSq.push_back( vRadius[i] * vRadius[i]);
2330
2331 std::vector<MathVector<dim> > vCircleCenter;
2332 for(size_t i = 0; i < vCenter.size(); ++i)
2333 vCircleCenter.push_back( *(vCenter[i]));
2334
2335
2336 // reset counter
2337 int numMarkedRefine = 0, numMarkedCoarse = 0;
2338
2339
2340 std::vector<MathVector<dim> > vCornerCoords;
2341 MathVector<dim> ElemCenter;
2342
2343 // loop elements for marking
2344 const_iterator iter = spF->template begin<TElem>();
2345 const_iterator iterEnd = spF->template end<TElem>();
2346 for(; iter != iterEnd; ++iter)
2347 {
2348 // get element
2349 TElem* elem = *iter;
2350
2351 std::vector<DoFIndex> ind;
2352 if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
2353 const number val = DoFRef(*spF, ind[0]) * DoFRef(*spPsi0, ind[0]);
2354
2355 // check for extra refinement
2356 bool bInCircle = false;
2357 if(!vRadiusSq.empty()){
2358 CollectCornerCoordinates(vCornerCoords, *elem, aaPos);
2359 AveragePositions<dim>(ElemCenter, vCornerCoords);
2360
2361 for(size_t i = 0; i < vRadiusSq.size(); ++i){
2362 if(VecDistanceSq(vCircleCenter[i], ElemCenter) <= vRadiusSq[i]){
2363 bInCircle = true;
2364 }
2365 }
2366 }
2367
2368 // marks for refinement
2369 if( val > minValueToRefine || bInCircle)
2370 {
2371 if(spF->dd()->multi_grid()->get_level(elem) < maxLevel)
2372 {
2373 refiner.mark(elem, RM_REFINE);
2374 numMarkedRefine++;
2375 }
2376
2377 continue;
2378 }
2379
2380 // marks for coarsening
2381 if( val < maxValueToCoarsen)
2382 {
2383 refiner.mark(elem, RM_COARSEN);
2384 numMarkedCoarse++;
2385 }
2386
2387 }
2388
2390 // print infos
2392
2393#ifdef UG_PARALLEL
2394 if(pcl::NumProcs() > 1)
2395 {
2396 UG_LOG(" >>> Marked on Proc "<<pcl::ProcRank()<<": refine: " << numMarkedRefine << ", coarsen : " << numMarkedCoarse << "\n");
2397 pcl::ProcessCommunicator com;
2398 int numMarkedRefineLocal = numMarkedRefine, numMarkedCoarseLocal = numMarkedCoarse;
2399 numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
2400 numMarkedCoarse = com.allreduce(numMarkedCoarseLocal, PCL_RO_SUM);
2401 }
2402#endif
2403
2404 UG_LOG(" >>> Marked refine: " << numMarkedRefine << ", coarsen: " << numMarkedCoarse << "\n" );
2405}
2406
2407
2408
2409template<typename TDomain>
2410void MarkDamage_OLD_AND_DEPRECATED( SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF,
2411 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spPsi0,
2412 IRefiner& refiner,
2413 number refineFrac, number coarseFrac,
2414 number avgRefineFactor, number avgCoarsenFactor,
2415 int maxLevel)
2416{
2417 PROFILE_FUNC_GROUP("Small Strain Mech");
2418
2419 static const int dim = TDomain::dim;
2420 typedef typename grid_dim_traits<dim>::element_type TElem;
2421 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
2422 const int fct = 0;
2423
2425 // statistic
2427
2428 number maxValue = std::numeric_limits<number>::min();
2429 number minValue = std::numeric_limits<number>::max();
2430 number sumValue = 0.0;
2431
2432 // loop all elems
2433 const size_t numElem = spF->size();
2434 for(size_t i = 0; i < numElem; ++i){
2435 const number val = (*spF)[i] * (*spPsi0)[i];
2436
2437 maxValue = std::max(maxValue, val);
2438 minValue = std::min(minValue, val);
2439 sumValue += val;
2440 }
2441 number avgValue = sumValue / numElem;
2442
2443
2444 UG_LOG(" +++ numElem: " << numElem << "\n");
2445 UG_LOG(" +++ maxValue: " << maxValue << "\n");
2446 UG_LOG(" +++ minValue: " << minValue << "\n");
2447 UG_LOG(" +++ avgValue: " << avgValue << "\n");
2448
2449
2451 // selection criteria
2453
2454 number minValueToRefine = std::min( refineFrac * maxValue, avgRefineFactor * avgValue);
2455 number maxValueToCoarsen = std::max( (1 + coarseFrac) * minValue, avgCoarsenFactor * avgValue);
2456
2457 UG_LOG(" +++ minValueToRefine: " << minValueToRefine << "\n");
2458 UG_LOG(" +++ maxValueToCoarsen: " << maxValueToCoarsen << "\n");
2459
2461 // marking
2463
2464 // reset counter
2465 int numMarkedRefine = 0, numMarkedCoarse = 0;
2466
2467 const_iterator iter = spF->template begin<TElem>();
2468 const_iterator iterEnd = spF->template end<TElem>();
2469
2470 // loop elements for marking
2471 for(; iter != iterEnd; ++iter)
2472 {
2473 // get element
2474 TElem* elem = *iter;
2475
2476 std::vector<DoFIndex> ind;
2477 if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
2478 const number val = DoFRef(*spF, ind[0]) * DoFRef(*spPsi0, ind[0]);
2479
2480 // marks for refinement
2481 if( val >= minValueToRefine)
2482 if(spF->dd()->multi_grid()->get_level(elem) < maxLevel)
2483 {
2484 refiner.mark(elem, RM_REFINE);
2485 numMarkedRefine++;
2486 }
2487
2488 // marks for coarsening
2489 if( val <= maxValueToCoarsen)
2490 {
2491 refiner.mark(elem, RM_COARSEN);
2492 numMarkedCoarse++;
2493 }
2494
2495 }
2496
2498 // print infos
2500
2501#ifdef UG_PARALLEL
2502 if(pcl::NumProcs() > 1)
2503 {
2504 UG_LOG(" +++ Marked for refinement on Proc "<<pcl::ProcRank()<<": " << numMarkedRefine << " Elements.\n");
2505 UG_LOG(" +++ Marked for coarsening on Proc "<<pcl::ProcRank()<<": " << numMarkedCoarse << " Elements.\n");
2506 pcl::ProcessCommunicator com;
2507 int numMarkedRefineLocal = numMarkedRefine, numMarkedCoarseLocal = numMarkedCoarse;
2508 numMarkedRefine = com.allreduce(numMarkedRefineLocal, PCL_RO_SUM);
2509 numMarkedCoarse = com.allreduce(numMarkedCoarseLocal, PCL_RO_SUM);
2510 }
2511#endif
2512
2513 UG_LOG(" +++ Marked for refinement: " << numMarkedRefine << " Elements.\n");
2514 UG_LOG(" +++ Marked for coarsening: " << numMarkedCoarse << " Elements.\n");
2515}
2516
2517
2518template<typename TDomain>
2519std::vector<number> DamageStatistic( SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF,
2520 SmartPtr<GridFunction<TDomain, CPUAlgebra> > spPsi0)
2521{
2522 PROFILE_FUNC_GROUP("Small Strain Mech");
2523
2524 static const int dim = TDomain::dim;
2525 typedef typename grid_dim_traits<dim>::element_type TElem;
2526 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
2527 const int fct = 0;
2528
2530 // statistic
2532
2533 number max_fPsi0 = std::numeric_limits<number>::min();
2534 number min_fPsi0 = std::numeric_limits<number>::max();
2535 number sum_fPsi0 = 0.0;
2536
2537 number max_l2_fPsi0 = std::numeric_limits<number>::min();
2538 number min_l2_fPsi0 = std::numeric_limits<number>::max();
2539 number sum_l2_fPsi0 = 0.0;
2540
2541 number max_l_fPsi0 = std::numeric_limits<number>::min();
2542 number min_l_fPsi0 = std::numeric_limits<number>::max();
2543 number sum_l_fPsi0 = 0.0;
2544
2545 // loop all elems
2546 const_iterator iter = spF->template begin<TElem>();
2547 const_iterator iterEnd = spF->template end<TElem>();
2548
2549 // loop elements for marking
2550 size_t numElem = 0;
2551 for(; iter != iterEnd; ++iter)
2552 {
2553 // get element
2554 TElem* elem = *iter; ++numElem;
2555
2556 std::vector<DoFIndex> ind;
2557 if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs");
2558
2559 const number f = DoFRef(*spF, ind[0]);
2560 const number psi0 = DoFRef(*spPsi0, ind[0]);
2561
2562 max_fPsi0 = std::max(max_fPsi0, f*psi0);
2563 min_fPsi0 = std::min(min_fPsi0, f*psi0);
2564 sum_fPsi0 += f*psi0;
2565
2566 const number l2 = ElementDiameterSq(*elem, *(spF->domain()) );
2567
2568 max_l2_fPsi0 = std::max(max_l2_fPsi0, l2*f*psi0);
2569 min_l2_fPsi0 = std::min(min_l2_fPsi0, l2*f*psi0);
2570 sum_l2_fPsi0 += l2*f*psi0;
2571
2572 const number l = std::sqrt(l2);
2573
2574 max_l_fPsi0 = std::max(max_l_fPsi0, l*f*psi0);
2575 min_l_fPsi0 = std::min(min_l_fPsi0, l*f*psi0);
2576 sum_l_fPsi0 += l*f*psi0;
2577 }
2578 number avg_fPsi0 = sum_fPsi0 / numElem;
2579 number avg_l2_fPsi0 = sum_l2_fPsi0 / numElem;
2580 number avg_l_fPsi0 = sum_l_fPsi0 / numElem;
2581
2582
2583 std::vector<number> vRes;
2584 vRes.push_back(max_fPsi0);
2585 vRes.push_back(min_fPsi0);
2586 vRes.push_back(avg_fPsi0);
2587
2588 vRes.push_back(max_l2_fPsi0);
2589 vRes.push_back(min_l2_fPsi0);
2590 vRes.push_back(avg_l2_fPsi0);
2591
2592 vRes.push_back(max_l_fPsi0);
2593 vRes.push_back(min_l_fPsi0);
2594 vRes.push_back(avg_l_fPsi0);
2595
2596 return vRes;
2597}
2598
2599
2600
2601template<typename TDomain>
2602void HadamardProd( SmartPtr<GridFunction<TDomain, CPUAlgebra> > spFPsi0,
2603 ConstSmartPtr<GridFunction<TDomain, CPUAlgebra> > spF,
2604 ConstSmartPtr<GridFunction<TDomain, CPUAlgebra> > spPsi0)
2605{
2606 if(spF->size() != spPsi0->size() || spF->size() != spFPsi0->size())
2607 UG_THROW("HadamardProd: Size mismatch");
2608
2609 for(size_t i = 0; i < spF->size(); ++i){
2610 (*spFPsi0)[i] = (*spF)[i] * (*spPsi0)[i];
2611 }
2612}
2613
2614
2615
2616
2617template<typename TDomain>
2618std::vector<number> MinMaxElementDiameter(SmartPtr<GridFunction<TDomain, CPUAlgebra> > spF)
2619{
2620 PROFILE_FUNC_GROUP("Small Strain Mech");
2621
2622 static const int dim = TDomain::dim;
2623 typedef typename grid_dim_traits<dim>::element_type TElem;
2624 typedef typename DoFDistribution::traits<TElem>::const_iterator const_iterator;
2625
2626 typename TDomain::grid_type& grid = *(spF->domain()->grid());
2627 typename TDomain::position_accessor_type& aaPos = spF->domain()->position_accessor();
2628
2629 // loop all elems
2630 const_iterator iter = spF->template begin<TElem>();
2631 const_iterator iterEnd = spF->template end<TElem>();
2632
2633 // loop elements for marking
2634 number max = 0.0;
2635 number min = std::numeric_limits<number>::max();
2636 for(; iter != iterEnd; ++iter){
2637 const number size = ElementDiameterSq(grid, aaPos, *iter);
2638 max = std::max(max, size);
2639 min = std::min(min, size);
2640 }
2641
2642#ifdef UG_PARALLEL
2643 // share value between all procs
2644 pcl::ProcessCommunicator com;
2645 max = com.allreduce(max, PCL_RO_MAX);
2646 min = com.allreduce(min, PCL_RO_MIN);
2647#endif
2648
2649 std::vector<number> vRes;
2650
2651 vRes.push_back( std::sqrt(min) );
2652 vRes.push_back( std::sqrt(max) );
2653
2654 return vRes;
2655
2656}
2657
2658
2659
2660
2661} // end namespace SmallStrainMechanics
2662}// namespace ug
2663
2664#endif /* __SMALL_STRAIN_MECHANICS__DAMAGE_IMPL_H_ */
p
parameterString p
Definition Biogas.lua:1
s
parameterString s
Definition Biogas.lua:2
pcl::ProcessCommunicator::allreduce
size_t allreduce(const size_t &t, pcl::ReduceOperation op) const
ug::IRefiner::mark
void mark(const TIterator &iterBegin, const TIterator &iterEnd, RefinementMark refMark=RM_REFINE)
ug::SmallStrainMechanics::DamageFunctionUpdater::write_debug
void write_debug(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spGF, std::string name, int call, int iter)
Definition damage_impl.h:1774
ug::SmallStrainMechanics::DamageFunctionUpdater::set_debug
void set_debug(SmartPtr< GridFunctionDebugWriter< TDomain, CPUAlgebra > > spDebugWriter)
Definition damage_impl.h:1767
ug::SmallStrainMechanics::DamageFunctionUpdater::solve
bool solve(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, SmartPtr< GridFunction< TDomain, CPUAlgebra > > spPsi0, const number beta, const number r, const number eps, const int maxIter, const number dampNewton)
Definition damage_impl.h:1597
ug::SmallStrainMechanics::DamageFunctionUpdater::set_disc_type
void set_disc_type(const std::string &type)
Definition damage_impl.h:1580
ug::SmallStrainMechanics::DamageFunctionUpdater::write_stencil_matrix_debug
void write_stencil_matrix_debug(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spGF, std::string name, int call)
Definition damage_impl.h:1794
ug::SmallStrainMechanics::RelativeDensityUpdater::write_debug
void write_debug(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spGF, std::string name, int call, int iter)
Definition damage_impl.h:2094
ug::SmallStrainMechanics::RelativeDensityUpdater::write_stencil_matrix_debug
void write_stencil_matrix_debug(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spGF, std::string name, int call)
Definition damage_impl.h:2114
ug::SmallStrainMechanics::RelativeDensityUpdater::solve
std::vector< number > solve(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spChi, SmartPtr< GridFunction< TDomain, CPUAlgebra > > spDrivingForce, const number betaStar, const number etaChiStar, const number chiMin, const number dt, const int p, const number rho_target, const number MassTol)
Definition damage_impl.h:1858
ug::SmallStrainMechanics::RelativeDensityUpdater::set_debug
void set_debug(SmartPtr< GridFunctionDebugWriter< TDomain, CPUAlgebra > > spDebugWriter)
Definition damage_impl.h:2087
ug::SmallStrainMechanics::RelativeDensityUpdater::set_disc_type
void set_disc_type(const std::string &type)
Definition damage_impl.h:1840
grid
SmartPtr< TGrid > grid()
ug::ElementDiameterSq
number ElementDiameterSq(const TElem &elem, TDomain &domain)
ug::CollectCornerCoordinates
void CollectCornerCoordinates(int base_object_id, std::vector< typename TDomain::position_type > &vCornerCoordsOut, GridObject &elem, const TDomain &domain, bool clearContainer)
ug::RM_COARSEN
RM_COARSEN
ug::RM_REFINE
RM_REFINE
PCL_RO_SUM
#define PCL_RO_SUM
pcl::ProcRank
int ProcRank()
PCL_RO_MAX
#define PCL_RO_MAX
PCL_RO_MIN
#define PCL_RO_MIN
pcl::NumProcs
int NumProcs()
ug::Invert
bool Invert(DenseMatrix< FixedArray2< double, 1, 1 > > &mat)
type
Variant::Type type()
ug::TrimString
UG_API std::string TrimString(const std::string &str)
UG_THROW
#define UG_THROW(msg)
UG_LOG
#define UG_LOG(msg)
number
double number
ug::ProjectPointToPlane
void ProjectPointToPlane(vector_t &vOut, const vector_t &v, const vector_t &p, const vector_t &n)
ug::VecTwoNormSq
vector_t::value_type VecTwoNormSq(const vector_t &v)
ug::VecTwoNorm
vector_t::value_type VecTwoNorm(const vector_t &v)
ug::VecDot
vector_t::value_type VecDot(const vector_t &v1, const vector_t &v2)
math_vector_functions.h
ug::SmallStrainMechanics::InitLaplacian_PartialIntegration
void InitLaplacian_PartialIntegration(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, std::vector< std::vector< number > > &vStencil, std::vector< std::vector< size_t > > &vIndex, int quadRuleType, bool fillElemSizeIntoVector=false)
Definition damage_impl.h:621
ug::SmallStrainMechanics::MarkDamage
void MarkDamage(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, SmartPtr< GridFunction< TDomain, CPUAlgebra > > spPsi0, IRefiner &refiner, number minValueToRefine, number maxValueToCoarsen, int maxLevel, const std::vector< MathVector< TDomain::dim, number > * > &vCenter, const std::vector< number > &vRadius)
Definition damage_impl.h:2293
ug::SmallStrainMechanics::MarkForAdaption_ValueRangeIndicator
void MarkForAdaption_ValueRangeIndicator(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spChi, IRefiner &refiner, number lowerValueToCoarsen, number minValueToRefine, number maxValueToRefine, number upperValueToCoarsen, number maxJumpDifferenceToNeighborsForCoarsen, number minJumpDifferenceToNeighborsForRefine, int maxLevel)
Definition damage_impl.h:2160
ug::SmallStrainMechanics::AveragePositions
void AveragePositions(MathVector< dim > &vCenter, const std::vector< MathVector< dim > > &vCornerCoords)
Definition damage_impl.h:44
ug::SmallStrainMechanics::MinMaxElementDiameter
std::vector< number > MinMaxElementDiameter(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF)
Definition damage_impl.h:2618
ug::SmallStrainMechanics::InitLaplacian_LeastSquares
void InitLaplacian_LeastSquares(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, std::vector< std::vector< number > > &vStencil, std::vector< std::vector< size_t > > &vIndex, bool fillElemSizeIntoVector=false)
Definition damage_impl.h:1242
ug::SmallStrainMechanics::HadamardProd
void HadamardProd(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spFPsi0, ConstSmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, ConstSmartPtr< GridFunction< TDomain, CPUAlgebra > > spPsi0)
Definition damage_impl.h:2602
ug::SmallStrainMechanics::CollectStencilNeighbors_NeumannZeroBND_IndexAndDistance
void CollectStencilNeighbors_NeumannZeroBND_IndexAndDistance(std::vector< typename grid_dim_traits< TDomain::dim >::element_type * > &vElem, std::vector< size_t > &vIndex, std::vector< MathVector< TDomain::dim > > &vDistance, typename grid_dim_traits< TDomain::dim >::element_type *elem, typename TDomain::grid_type &grid, typename TDomain::position_accessor_type &aaPos, SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, bool fillElemSizeIntoVector=false)
Definition damage_impl.h:297
ug::SmallStrainMechanics::CollectSurfaceNeighbors
void CollectSurfaceNeighbors(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, typename grid_dim_traits< TDomain::dim >::element_type *elem, std::vector< typename grid_dim_traits< TDomain::dim >::element_type * > &vNeighbors)
Definition damage_impl.h:62
ug::SmallStrainMechanics::InitLaplacian_TaylorExpansion
void InitLaplacian_TaylorExpansion(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, std::vector< std::vector< number > > &vStencil, std::vector< std::vector< size_t > > &vIndex, bool fillElemSizeIntoVector=false)
Definition damage_impl.h:1035
ug::SmallStrainMechanics::MarkDamage_OLD_AND_DEPRECATED
void MarkDamage_OLD_AND_DEPRECATED(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, SmartPtr< GridFunction< TDomain, CPUAlgebra > > spPsi0, IRefiner &refiner, number refineFrac, number coarseFrac, number avgRefineFactor, number avgCoarsenFactor, int maxLevel)
Definition damage_impl.h:2410
ug::SmallStrainMechanics::InitLaplacian_TaylorDirect
void InitLaplacian_TaylorDirect(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, std::vector< std::vector< number > > &vStencil, std::vector< std::vector< size_t > > &vIndex, bool fillElemSizeIntoVector=false)
Definition damage_impl.h:1417
ug::SmallStrainMechanics::DamageStatistic
std::vector< number > DamageStatistic(SmartPtr< GridFunction< TDomain, CPUAlgebra > > spF, SmartPtr< GridFunction< TDomain, CPUAlgebra > > spPsi0)
Definition damage_impl.h:2519
ug::DoFRef
const number & DoFRef(const TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
ug::VecDistanceSq
TVector::value_type VecDistanceSq(const TVector &v1, const TVector &v2, const TMatrix &M)
ug::VecScaleAdd
void VecScaleAdd(double &dest, double alpha1, const double &v1, double alpha2, const double &v2)
PROFILE_FUNC_GROUP
#define PROFILE_FUNC_GROUP(groups)
PROFILE_BEGIN_GROUP
#define PROFILE_BEGIN_GROUP(name, groups)
ug::SmallStrainMechanics::contrained_dim_traits
Definition damage.h:50
ug::grid_dim_traits