Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

# This file is part of the pyMOR project (http://www.pymor.org). 

# Copyright Holders: Rene Milk, Stephan Rave, Felix Schindler 

# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) 

 

from __future__ import absolute_import, division, print_function 

 

import numpy as np 

 

from pymor.core import abstractmethod, abstractstaticmethod 

from pymor.core.cache import CacheableInterface, cached 

from pymor.domaindescriptions import BoundaryType 

from pymor.grids.defaultimpl import (ConformalTopologicalGridDefaultImplementations, 

                                     SimpleReferenceElementDefaultImplementations, 

                                     AffineGridDefaultImplementations,) 

 

 

class ConformalTopologicalGridInterface(ConformalTopologicalGridDefaultImplementations, CacheableInterface): 

    '''A conformal topological grid. 

 

    The grid is completely determined via the subentity relation given by 

    :meth:`~ConformalTopologicalGridInterface.subentities`. In addition, 

    only :meth:`~ConformalTopologicalGridInterface.size` has to be 

    implemented, cached default implementations for all other methods are 

    provided by :class:`~pymor.grids.defaultimpl.ConformalTopologicalGridDefaultImplementations`. 

 

    Attributes 

    ---------- 

    dim 

        The dimension of the grid. 

    ''' 

 

    @abstractmethod 

    def size(self, codim): 

        '''The number of entities of codimension `codim`.''' 

        pass 

 

    @abstractmethod 

    def subentities(self, codim, subentity_codim): 

        '''`retval[e,s]` is the global index of the `s`-th codim-`subentity_codim` 

        subentity of the codim-`codim` entity with global index `e`. 

 

        Only `subentities(codim, codim+1)` has to be implemented; a default 

        implementation is provided which evaluates 

        `subentities(codim, subentity_codim)` by computing the 

        transitive closure of `subentities(codim, codim+1)`. 

        ''' 

        return self._subentities(codim, subentity_codim) 

 

    def superentities(self, codim, superentity_codim): 

        '''`retval[e,s]` is the global index of the `s`-th codim-`superentity_codim` 

        superentity of the codim-`codim` entity with global index `e`. 

 

        `retval[e]` is sorted by global index. 

 

        The default implementation is to compute the result from 

        `subentities(superentity_codim, codim)`. 

        ''' 

        return self._superentities(codim, superentity_codim) 

 

    def superentity_indices(self, codim, superentity_codim): 

        '''`retval[e,s]` is the local index of the codim-`codim` entity `e` 

        in the codim-`superentity_codim` superentity `superentities(codim, superentity_codim)[e,s].` 

        ''' 

        return self._superentity_indices(codim, superentity_codim) 

 

    def neighbours(self, codim, neighbour_codim, intersection_codim=None): 

        '''`retval[e,n]` is the global index of the `n`-th codim-`neighbour_codim` entitiy of the 

        codim-`codim` entity `e` that shares with `e` a subentity of codimension `intersection_codim`. 

 

        If `intersection_codim == None`, it is set to `codim + 1` if `codim == neighbour_codim` 

        and to `min(codim, neighbour_codim)` otherwise. 

 

        The default implementation is to compute the result from 

        `subentities(codim, intersection_codim)` and 

        `superentities(intersection_codim, neihbour_codim)`. 

        ''' 

        return self._neighbours(codim, neighbour_codim, intersection_codim) 

 

    def boundary_mask(self, codim): 

        '''`retval[e]` is true iff the codim-`codim` entity with global index 

        `e` is a boundary entity. 

 

        By definition, a codim-1 entity is a boundary entity if it has only one 

        codim-0 superentity. For `codim != 1`, a codim-`codim` entity is a 

        boundary entity if it has a codim-1 sub/super-entity. 

        ''' 

        return self._boundary_mask(codim) 

 

    def boundaries(self, codim): 

        '''Returns the global indices of all codim-`codim` boundary entities. 

 

        By definition, a codim-1 entity is a boundary entity if it has only one 

        codim-0 superentity. For `codim != 1`, a codim-`codim` entity is a 

        boundary entity if it has a codim-1 sub/super-entity. 

        ''' 

        return self._boundaries(codim) 

 

    @abstractstaticmethod 

    def test_instances(): 

        '''Returns a list of Grid instances suitable to be run through our test cases.''' 

        pass 

 

 

class ReferenceElementInterface(SimpleReferenceElementDefaultImplementations, CacheableInterface): 

    '''Defines a reference element with the property that each of its subentities is of the same type. 

 

    I.e. a three-dimensional reference element cannot have triangles and rectangles as faces at the 

    same time. 

 

    Attributes 

    ---------- 

    dim 

        The dimension of the reference element 

    volume 

        The volume of the reference element 

    ''' 

 

    dim = None 

    volume = None 

 

    @abstractmethod 

    def size(self, codim): 

        'Number of subentites of codimension `codim`.' 

 

    @abstractmethod 

    def subentities(self, codim, subentity_codim): 

        '''`subentities(c,sc)[i,j]` is, with respect to the indexing inside the 

        reference element, the index of the `j`-th codim-`subentity_codim` 

        subentity of the `i`-th codim-`codim` subentity of the reference element. 

        ''' 

        pass 

 

    @abstractmethod 

    def subentity_embedding(self, subentity_codim): 

        '''Returns a tuple `(A, B)` which defines the embedding of the codim-`subentity_codim` 

        subentities into the reference element. 

 

        For `subentity_codim > 1', the embedding is by default given recursively via 

        `subentity_embedding(subentity_codim - 1)` and 

        `sub_reference_element(subentity_codim - 1).subentity_embedding(1)` choosing always 

        the superentity with smallest index. 

        ''' 

        return self._subentity_embedding(subentity_codim) 

 

    @abstractmethod 

    def sub_reference_element(self, codim): 

        '''Returns the reference relement of the codim-`codim` subentities.''' 

        return self._sub_reference_element(codim) 

 

    def __call__(self, codim): 

        '''Returns the reference relement of the codim-`codim` subentities.''' 

        return self.sub_reference_element(codim) 

 

    @abstractmethod 

    def unit_outer_normals(self): 

        '''`retval[e]` is the unit outer-normal vector to the codim-1 subentity 

        with index `e`. 

        ''' 

        pass 

 

    @abstractmethod 

    def center(self): 

        '''Coordinates of the barycenter.''' 

        pass 

 

    @abstractmethod 

    def mapped_diameter(self, A): 

        '''The diameter of the reference element after tranforming it with the 

        matrix `A` (vectorized). 

        ''' 

        pass 

 

    @abstractmethod 

    def quadrature(self, order=None, npoints=None, quadrature_type='default'): 

        '''Returns tuple `(P, W)` where `P` is an array of quadrature points with 

        corresponding weights `W`. 

 

        The quadrature is of order `order` or has `npoints` integration points. 

        ''' 

        pass 

 

    @abstractmethod 

    def quadrature_info(self): 

        '''Returns a tuple of dicts `(O, N)` where `O[quadrature_type]` is a list 

        of orders which are implemented for `quadrature_type` and `N[quadrature_type]` 

        is a list of the corrsponding numbers of integration points. 

        ''' 

        pass 

 

    def quadrature_types(self): 

        o, _ = self.quadrature_info() 

        return frozenset(o.keys()) 

 

 

class AffineGridInterface(AffineGridDefaultImplementations, ConformalTopologicalGridInterface): 

    '''Topological grid with geometry where each codim-0 entity is affinely mapped to the same |ReferenceElement|. 

 

    The grid is completely determined via the subentity relation given by 

    :meth:`~AffineGridInterface.subentities` and the embeddings given by 

    :meth:`~AffineGridInterface.embeddings`.  In addition, 

    only :meth:`~ConformalTopologicalGridInterface.size` and :meth:`~AffineGridInterface.reference_element` 

    have to be implemented. Cached default implementations for all other methods are 

    provided by :class:`~pymor.grids.defaultimpl.AffineGridDefaultImplementations`. 

 

    Attributes 

    ---------- 

    dim_outer 

        The dimension of the space into which the grid is embedded. 

    ''' 

 

    dim_outer = None 

 

    @abstractmethod 

    def reference_element(self, codim): 

        '''The |ReferenceElement| of the codim-`codim` entities.''' 

        pass 

 

    @abstractmethod 

    def subentities(self, codim, subentity_codim): 

        '''`retval[e,s]` is the global index of the `s`-th codim-`subentity_codim` 

        subentity of the codim-`codim` entity with global index `e`. 

 

        The ordering of `subentities(0, subentity_codim)[e]` has to correspond, w.r.t. 

        the embedding of `e`, to the local ordering inside the reference element. 

 

        For `codim > 0`, we provide a default implementation by calculating the 

        subentites of `e` as follows: 

 

            1. Find the `codim-1` parent entity `e_0` of `e` with minimal global index 

            2. Lookup the local indicies of the subentites of `e` inside `e_0` using the reference element. 

            3. Map these local indicies to global indicies using `subentities(codim - 1, subentity_codim)`. 

 

        This procedures assures that `subentities(codim, subentity_codim)[e]` 

        has the right ordering w.r.t. the embedding determined by `e_0`, which 

        agrees with what is returned by `embeddings(codim)` 

        ''' 

        return self._subentities(codim, subentity_codim) 

 

    @abstractmethod 

    def embeddings(self, codim): 

        '''Returns tuple `(A, B)` where `A[e]` and `B[e]` are the linear part 

        and the translation part of the map from the reference element of `e` 

        to `e`. 

 

        For `codim > 0`, we provide a default implementation by 

        taking the embedding of the codim-1 parent entity `e_0` of `e` with 

        lowest global index and composing it with the subentity_embedding of `e` 

        into `e_0` determined by the reference element. 

        ''' 

        return self._embeddings(codim) 

 

    def jacobian_inverse_transposed(self, codim): 

        '''`retval[e]` is the transposed (pseudo-)inverse of the jacobian of `embeddings(codim)[e]`. 

        ''' 

        return self._jacobian_inverse_transposed(codim) 

 

    def integration_elements(self, codim): 

        '''`retval[e]` is given as `sqrt(det(A^T*A))`, where `A = embeddings(codim)[0][e]`.''' 

        return self._integration_elements(codim) 

 

    def volumes(self, codim): 

        '''`retval[e]` is the (dim-`codim`)-dimensional volume of the codim-`codim` entity with global index `e`.''' 

        return self._volumes(codim) 

 

    def volumes_inverse(self, codim): 

        '''`retval[e] = 1 / volumes(codim)[e]`.''' 

        return self._volumes_inverse(codim) 

 

    def unit_outer_normals(self): 

        '''`retval[e,i]` is the unit outer normal to the i-th codim-1 subentity 

        of the codim-0 entitiy with global index `e`. 

        ''' 

        return self._unit_outer_normals() 

 

    def centers(self, codim): 

        '''`retval[e]` is the barycenter of the codim-`codim` entity with global index `e`.''' 

        return self._centers(codim) 

 

    def diameters(self, codim): 

        '''`retval[e]` is the diameter of the codim-`codim` entity with global index `e`.''' 

        return self._diameters(codim) 

 

    def quadrature_points(self, codim, order=None, npoints=None, quadrature_type='default'): 

        '''`retval[e]` is an array of quadrature points in global coordinates 

        for the codim-`codim` entity with global index `e`. 

 

        The quadrature is of order `order` or has `npoints` integration points. To 

        integrate a function `f` over `e` one has to form :: 

 

            np.dot(f(quadrature_points(codim, order)[e]), reference_element(codim).quadrature(order)[1]) * integration_elements(codim)[e].  # NOQA 

        ''' 

        return self._quadrature_points(codim, order, npoints, quadrature_type) 

 

 

class BoundaryInfoInterface(CacheableInterface): 

    '''Provides |BoundaryTypes| for the boundaries of a given |ConformalTopologicalGrid|. 

 

    For every |BoundaryType| and codimension a mask is provided, marking grid entities 

    of the respective type and codimension by their global index. 

 

    Attributes 

    ---------- 

    boundary_types 

        set of all |BoundaryTypes| the grid has. 

    ''' 

 

    boundary_types = frozenset() 

 

    def mask(self, boundary_type, codim): 

        '''retval[i] is `True` if the codim-`codim` entity of global index `i` is 

        associated to the |BoundaryType| `boundary_type`. 

        ''' 

        raise ValueError('Has no boundary_type "{}"'.format(boundary_type)) 

 

    def unique_boundary_type_mask(self, codim): 

        '''retval[i] is `True` if the codim-`codim` entity of global index `i` is 

        associated to one and only one |BoundaryType|. 

        ''' 

        return np.less_equal(sum(self.mask(bt, codim=codim).astype(np.int) for bt in self.boundary_types), 1) 

 

    def no_boundary_type_mask(self, codim): 

        '''retval[i] is `True` if the codim-`codim` entity of global index `i` is 

        associated to no |BoundaryType|. 

        ''' 

        return np.equal(sum(self.mask(bt, codim=codim).astype(np.int) for bt in self.boundary_types), 0) 

 

    def check_boundary_types(self, assert_unique_type=[1], assert_some_type=[]): 

        if assert_unique_type: 

            for codim in assert_unique_type: 

                assert np.all(self.unique_boundary_type_mask(codim)) 

        if assert_some_type: 

            for codim in assert_some_type: 

                assert not np.any(self.no_boundary_type_mask(codim)) 

 

    @property 

    def has_dirichlet(self): 

        return BoundaryType('dirichlet') in self.boundary_types 

 

    @property 

    def has_neumann(self): 

        return BoundaryType('neumann') in self.boundary_types 

 

    @property 

    def has_only_dirichlet(self): 

        return self.boundary_types == {BoundaryType('dirichlet')} 

 

    @property 

    def has_only_neumann(self): 

        return self.boundary_types == {BoundaryType('neumann')} 

 

    @property 

    def has_only_dirichletneumann(self): 

        return self.boundary_types <= {BoundaryType('dirichlet'), BoundaryType('neumann')} 

 

    def dirichlet_mask(self, codim): 

        return self.mask(BoundaryType('dirichlet'), codim) 

 

    def neumann_mask(self, codim): 

        return self.mask(BoundaryType('neumann'), codim) 

 

    @cached 

    def _dirichlet_boundaries(self, codim): 

        return np.where(self.dirichlet_mask(codim))[0].astype('int32') 

 

    def dirichlet_boundaries(self, codim): 

        return self._dirichlet_boundaries(codim) 

 

    @cached 

    def _neumann_boundaries(self, codim): 

        return np.where(self.neumann_mask(codim))[0].astype('int32') 

 

    def neumann_boundaries(self, codim): 

        return self._neumann_boundaries(codim)