1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petscdmplex.h> 3 #include <petscblaslapack.h> 4 5 PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]); 6 7 struct _n_Petsc1DNodeFamily 8 { 9 PetscInt refct; 10 PetscDTNodeType nodeFamily; 11 PetscReal gaussJacobiExp; 12 PetscInt nComputed; 13 PetscReal **nodesets; 14 PetscBool endpoints; 15 }; 16 17 /* users set node families for PETSCDUALSPACELAGRANGE with just the inputs to this function, but internally we create 18 * an object that can cache the computations across multiple dual spaces */ 19 static PetscErrorCode Petsc1DNodeFamilyCreate(PetscDTNodeType family, PetscReal gaussJacobiExp, PetscBool endpoints, Petsc1DNodeFamily *nf) 20 { 21 Petsc1DNodeFamily f; 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 ierr = PetscNew(&f);CHKERRQ(ierr); 26 switch (family) { 27 case PETSCDTNODES_GAUSSJACOBI: 28 case PETSCDTNODES_EQUISPACED: 29 f->nodeFamily = family; 30 break; 31 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family"); 32 } 33 f->endpoints = endpoints; 34 f->gaussJacobiExp = 0.; 35 if (family == PETSCDTNODES_GAUSSJACOBI) { 36 if (gaussJacobiExp <= -1.) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Gauss-Jacobi exponent must be > -1.\n"); 37 f->gaussJacobiExp = gaussJacobiExp; 38 } 39 f->refct = 1; 40 *nf = f; 41 PetscFunctionReturn(0); 42 } 43 44 static PetscErrorCode Petsc1DNodeFamilyReference(Petsc1DNodeFamily nf) 45 { 46 PetscFunctionBegin; 47 if (nf) nf->refct++; 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode Petsc1DNodeFamilyDestroy(Petsc1DNodeFamily *nf) 52 { 53 PetscInt i, nc; 54 PetscErrorCode ierr; 55 56 PetscFunctionBegin; 57 if (!(*nf)) PetscFunctionReturn(0); 58 if (--(*nf)->refct > 0) { 59 *nf = NULL; 60 PetscFunctionReturn(0); 61 } 62 nc = (*nf)->nComputed; 63 for (i = 0; i < nc; i++) { 64 ierr = PetscFree((*nf)->nodesets[i]);CHKERRQ(ierr); 65 } 66 ierr = PetscFree((*nf)->nodesets);CHKERRQ(ierr); 67 ierr = PetscFree(*nf);CHKERRQ(ierr); 68 *nf = NULL; 69 PetscFunctionReturn(0); 70 } 71 72 static PetscErrorCode Petsc1DNodeFamilyGetNodeSets(Petsc1DNodeFamily f, PetscInt degree, PetscReal ***nodesets) 73 { 74 PetscInt nc; 75 PetscErrorCode ierr; 76 77 PetscFunctionBegin; 78 nc = f->nComputed; 79 if (degree >= nc) { 80 PetscInt i, j; 81 PetscReal **new_nodesets; 82 PetscReal *w; 83 84 ierr = PetscMalloc1(degree + 1, &new_nodesets);CHKERRQ(ierr); 85 ierr = PetscArraycpy(new_nodesets, f->nodesets, nc);CHKERRQ(ierr); 86 ierr = PetscFree(f->nodesets);CHKERRQ(ierr); 87 f->nodesets = new_nodesets; 88 ierr = PetscMalloc1(degree + 1, &w);CHKERRQ(ierr); 89 for (i = nc; i < degree + 1; i++) { 90 ierr = PetscMalloc1(i + 1, &(f->nodesets[i]));CHKERRQ(ierr); 91 if (!i) { 92 f->nodesets[i][0] = 0.5; 93 } else { 94 switch (f->nodeFamily) { 95 case PETSCDTNODES_EQUISPACED: 96 if (f->endpoints) { 97 for (j = 0; j <= i; j++) f->nodesets[i][j] = (PetscReal) j / (PetscReal) i; 98 } else { 99 /* these nodes are at the centroids of the small simplices created by the equispaced nodes that include 100 * the endpoints */ 101 for (j = 0; j <= i; j++) f->nodesets[i][j] = ((PetscReal) j + 0.5) / ((PetscReal) i + 1.); 102 } 103 break; 104 case PETSCDTNODES_GAUSSJACOBI: 105 if (f->endpoints) { 106 ierr = PetscDTGaussLobattoJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w);CHKERRQ(ierr); 107 } else { 108 ierr = PetscDTGaussJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w);CHKERRQ(ierr); 109 } 110 break; 111 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family"); 112 } 113 } 114 } 115 ierr = PetscFree(w);CHKERRQ(ierr); 116 f->nComputed = degree + 1; 117 } 118 *nodesets = f->nodesets; 119 PetscFunctionReturn(0); 120 } 121 122 /* http://arxiv.org/abs/2002.09421 for details */ 123 static PetscErrorCode PetscNodeRecursive_Internal(PetscInt dim, PetscInt degree, PetscReal **nodesets, PetscInt tup[], PetscReal node[]) 124 { 125 PetscReal w; 126 PetscInt i, j; 127 PetscErrorCode ierr; 128 129 PetscFunctionBeginHot; 130 w = 0.; 131 if (dim == 1) { 132 node[0] = nodesets[degree][tup[0]]; 133 node[1] = nodesets[degree][tup[1]]; 134 } else { 135 for (i = 0; i < dim + 1; i++) node[i] = 0.; 136 for (i = 0; i < dim + 1; i++) { 137 PetscReal wi = nodesets[degree][degree-tup[i]]; 138 139 for (j = 0; j < dim+1; j++) tup[dim+1+j] = tup[j+(j>=i)]; 140 ierr = PetscNodeRecursive_Internal(dim-1,degree-tup[i],nodesets,&tup[dim+1],&node[dim+1]);CHKERRQ(ierr); 141 for (j = 0; j < dim+1; j++) node[j+(j>=i)] += wi * node[dim+1+j]; 142 w += wi; 143 } 144 for (i = 0; i < dim+1; i++) node[i] /= w; 145 } 146 PetscFunctionReturn(0); 147 } 148 149 /* compute simplex nodes for the biunit simplex from the 1D node family */ 150 static PetscErrorCode Petsc1DNodeFamilyComputeSimplexNodes(Petsc1DNodeFamily f, PetscInt dim, PetscInt degree, PetscReal points[]) 151 { 152 PetscInt *tup; 153 PetscInt k; 154 PetscInt npoints; 155 PetscReal **nodesets = NULL; 156 PetscInt worksize; 157 PetscReal *nodework; 158 PetscInt *tupwork; 159 PetscErrorCode ierr; 160 161 PetscFunctionBegin; 162 if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have non-negative dimension\n"); 163 if (degree < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have non-negative degree\n"); 164 if (!dim) PetscFunctionReturn(0); 165 ierr = PetscCalloc1(dim+2, &tup);CHKERRQ(ierr); 166 k = 0; 167 ierr = PetscDTBinomialInt(degree + dim, dim, &npoints);CHKERRQ(ierr); 168 ierr = Petsc1DNodeFamilyGetNodeSets(f, degree, &nodesets);CHKERRQ(ierr); 169 worksize = ((dim + 2) * (dim + 3)) / 2; 170 ierr = PetscMalloc2(worksize, &nodework, worksize, &tupwork);CHKERRQ(ierr); 171 /* loop over the tuples of length dim with sum at most degree */ 172 for (k = 0; k < npoints; k++) { 173 PetscInt i; 174 175 /* turn thm into tuples of length dim + 1 with sum equal to degree (barycentric indice) */ 176 tup[0] = degree; 177 for (i = 0; i < dim; i++) { 178 tup[0] -= tup[i+1]; 179 } 180 switch(f->nodeFamily) { 181 case PETSCDTNODES_EQUISPACED: 182 /* compute equispaces nodes on the unit reference triangle */ 183 if (f->endpoints) { 184 for (i = 0; i < dim; i++) { 185 points[dim*k + i] = (PetscReal) tup[i+1] / (PetscReal) degree; 186 } 187 } else { 188 for (i = 0; i < dim; i++) { 189 /* these nodes are at the centroids of the small simplices created by the equispaced nodes that include 190 * the endpoints */ 191 points[dim*k + i] = ((PetscReal) tup[i+1] + 1./(dim+1.)) / (PetscReal) (degree + 1.); 192 } 193 } 194 break; 195 default: 196 /* compute equispaces nodes on the barycentric reference triangle (the trace on the first dim dimensions are the 197 * unit reference triangle nodes */ 198 for (i = 0; i < dim + 1; i++) tupwork[i] = tup[i]; 199 ierr = PetscNodeRecursive_Internal(dim, degree, nodesets, tupwork, nodework);CHKERRQ(ierr); 200 for (i = 0; i < dim; i++) points[dim*k + i] = nodework[i + 1]; 201 break; 202 } 203 ierr = PetscDualSpaceLatticePointLexicographic_Internal(dim, degree, &tup[1]);CHKERRQ(ierr); 204 } 205 /* map from unit simplex to biunit simplex */ 206 for (k = 0; k < npoints * dim; k++) points[k] = points[k] * 2. - 1.; 207 ierr = PetscFree2(nodework, tupwork);CHKERRQ(ierr); 208 ierr = PetscFree(tup);CHKERRQ(ierr); 209 PetscFunctionReturn(0); 210 } 211 212 /* If we need to get the dofs from a mesh point, or add values into dofs at a mesh point, and there is more than one dof 213 * on that mesh point, we have to be careful about getting/adding everything in the right place. 214 * 215 * With nodal dofs like PETSCDUALSPACELAGRANGE makes, the general approach to calculate the value of dofs associate 216 * with a node A is 217 * - transform the node locations x(A) by the map that takes the mesh point to its reorientation, x' = phi(x(A)) 218 * - figure out which node was originally at the location of the transformed point, A' = idx(x') 219 * - if the dofs are not scalars, figure out how to represent the transformed dofs in terms of the basis 220 * of dofs at A' (using pushforward/pullback rules) 221 * 222 * The one sticky point with this approach is the "A' = idx(x')" step: trying to go from real valued coordinates 223 * back to indices. I don't want to rely on floating point tolerances. Additionally, PETSCDUALSPACELAGRANGE may 224 * eventually support quasi-Lagrangian dofs, which could involve quadrature at multiple points, so the location "x(A)" 225 * would be ambiguous. 226 * 227 * So each dof gets an integer value coordinate (nodeIdx in the structure below). The choice of integer coordinates 228 * is somewhat arbitrary, as long as all of the relevant symmetries of the mesh point correspond to *permutations* of 229 * the integer coordinates, which do not depend on numerical precision. 230 * 231 * So 232 * 233 * - DMPlexGetTransitiveClosure_Internal() tells me how an orientation turns into a permutation of the vertices of a 234 * mesh point 235 * - The permutation of the vertices, and the nodeIdx values assigned to them, tells what permutation in index space 236 * is associated with the orientation 237 * - I uses that permutation to get xi' = phi(xi(A)), the integer coordinate of the transformed dof 238 * - I can without numerical issues compute A' = idx(xi') 239 * 240 * Here are some examples of how the process works 241 * 242 * - With a triangle: 243 * 244 * The triangle has the following integer coordinates for vertices, taken from the barycentric triangle 245 * 246 * closure order 2 247 * nodeIdx (0,0,1) 248 * \ 249 * + 250 * |\ 251 * | \ 252 * | \ 253 * | \ closure order 1 254 * | \ / nodeIdx (0,1,0) 255 * +-----+ 256 * \ 257 * closure order 0 258 * nodeIdx (1,0,0) 259 * 260 * If I do DMPlexGetTransitiveClosure_Internal() with orientation 1, the vertices would appear 261 * in the order (1, 2, 0) 262 * 263 * If I list the nodeIdx of each vertex in closure order for orientation 0 (0, 1, 2) and orientation 1 (1, 2, 0), I 264 * see 265 * 266 * orientation 0 | orientation 1 267 * 268 * [0] (1,0,0) [1] (0,1,0) 269 * [1] (0,1,0) [2] (0,0,1) 270 * [2] (0,0,1) [0] (1,0,0) 271 * A B 272 * 273 * In other words, B is the result of a row permutation of A. But, there is also 274 * a column permutation that accomplishes the same result, (2,0,1). 275 * 276 * So if a dof has nodeIdx coordinate (a,b,c), after the transformation its nodeIdx coordinate 277 * is (c,a,b), and the transformed degree of freedom will be a linear combination of dofs 278 * that originally had coordinate (c,a,b). 279 * 280 * - With a quadrilateral: 281 * 282 * The quadrilateral has the following integer coordinates for vertices, taken from concatenating barycentric 283 * coordinates for two segments: 284 * 285 * closure order 3 closure order 2 286 * nodeIdx (1,0,0,1) nodeIdx (0,1,0,1) 287 * \ / 288 * +----+ 289 * | | 290 * | | 291 * +----+ 292 * / \ 293 * closure order 0 closure order 1 294 * nodeIdx (1,0,1,0) nodeIdx (0,1,1,0) 295 * 296 * If I do DMPlexGetTransitiveClosure_Internal() with orientation 1, the vertices would appear 297 * in the order (1, 2, 3, 0) 298 * 299 * If I list the nodeIdx of each vertex in closure order for orientation 0 (0, 1, 2, 3) and 300 * orientation 1 (1, 2, 3, 0), I see 301 * 302 * orientation 0 | orientation 1 303 * 304 * [0] (1,0,1,0) [1] (0,1,1,0) 305 * [1] (0,1,1,0) [2] (0,1,0,1) 306 * [2] (0,1,0,1) [3] (1,0,0,1) 307 * [3] (1,0,0,1) [0] (1,0,1,0) 308 * A B 309 * 310 * The column permutation that accomplishes the same result is (3,2,0,1). 311 * 312 * So if a dof has nodeIdx coordinate (a,b,c,d), after the transformation its nodeIdx coordinate 313 * is (d,c,a,b), and the transformed degree of freedom will be a linear combination of dofs 314 * that originally had coordinate (d,c,a,b). 315 * 316 * Previously PETSCDUALSPACELAGRANGE had hardcoded symmetries for the triangle and quadrilateral, 317 * but this approach will work for any polytope, such as the wedge (triangular prism). 318 */ 319 struct _n_PetscLagNodeIndices 320 { 321 PetscInt refct; 322 PetscInt nodeIdxDim; 323 PetscInt nodeVecDim; 324 PetscInt nNodes; 325 PetscInt *nodeIdx; /* for each node an index of size nodeIdxDim */ 326 PetscReal *nodeVec; /* for each node a vector of size nodeVecDim */ 327 PetscInt *perm; /* if these are vertices, perm takes DMPlex point index to closure order; 328 if these are nodes, perm lists nodes in index revlex order */ 329 }; 330 331 /* this is just here so I can access the values in tests/ex1.c outside the library */ 332 PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices ni, PetscInt *nodeIdxDim, PetscInt *nodeVecDim, PetscInt *nNodes, const PetscInt *nodeIdx[], const PetscReal *nodeVec[]) 333 { 334 PetscFunctionBegin; 335 *nodeIdxDim = ni->nodeIdxDim; 336 *nodeVecDim = ni->nodeVecDim; 337 *nNodes = ni->nNodes; 338 *nodeIdx = ni->nodeIdx; 339 *nodeVec = ni->nodeVec; 340 PetscFunctionReturn(0); 341 } 342 343 static PetscErrorCode PetscLagNodeIndicesReference(PetscLagNodeIndices ni) 344 { 345 PetscFunctionBegin; 346 if (ni) ni->refct++; 347 PetscFunctionReturn(0); 348 } 349 350 static PetscErrorCode PetscLagNodeIndicesDuplicate(PetscLagNodeIndices ni, PetscLagNodeIndices *niNew) 351 { 352 PetscErrorCode ierr; 353 354 PetscFunctionBegin; 355 ierr = PetscNew(niNew);CHKERRQ(ierr); 356 (*niNew)->refct = 1; 357 (*niNew)->nodeIdxDim = ni->nodeIdxDim; 358 (*niNew)->nodeVecDim = ni->nodeVecDim; 359 (*niNew)->nNodes = ni->nNodes; 360 ierr = PetscMalloc1(ni->nNodes * ni->nodeIdxDim, &((*niNew)->nodeIdx));CHKERRQ(ierr); 361 ierr = PetscArraycpy((*niNew)->nodeIdx, ni->nodeIdx, ni->nNodes * ni->nodeIdxDim);CHKERRQ(ierr); 362 ierr = PetscMalloc1(ni->nNodes * ni->nodeVecDim, &((*niNew)->nodeVec));CHKERRQ(ierr); 363 ierr = PetscArraycpy((*niNew)->nodeVec, ni->nodeVec, ni->nNodes * ni->nodeVecDim);CHKERRQ(ierr); 364 (*niNew)->perm = NULL; 365 PetscFunctionReturn(0); 366 } 367 368 static PetscErrorCode PetscLagNodeIndicesDestroy(PetscLagNodeIndices *ni) 369 { 370 PetscErrorCode ierr; 371 372 PetscFunctionBegin; 373 if (!(*ni)) PetscFunctionReturn(0); 374 if (--(*ni)->refct > 0) { 375 *ni = NULL; 376 PetscFunctionReturn(0); 377 } 378 ierr = PetscFree((*ni)->nodeIdx);CHKERRQ(ierr); 379 ierr = PetscFree((*ni)->nodeVec);CHKERRQ(ierr); 380 ierr = PetscFree((*ni)->perm);CHKERRQ(ierr); 381 ierr = PetscFree(*ni);CHKERRQ(ierr); 382 *ni = NULL; 383 PetscFunctionReturn(0); 384 } 385 386 /* The vertices are given nodeIdx coordinates (e.g. the corners of the barycentric triangle). Those coordinates are 387 * in some other order, and to understand the effect of different symmetries, we need them to be in closure order. 388 * 389 * If sortIdx is PETSC_FALSE, the coordinates are already in revlex order, otherwise we must sort them 390 * to that order before we do the real work of this function, which is 391 * 392 * - mark the vertices in closure order 393 * - sort them in revlex order 394 * - use the resulting permutation to list the vertex coordinates in closure order 395 */ 396 static PetscErrorCode PetscLagNodeIndicesComputeVertexOrder(DM dm, PetscLagNodeIndices ni, PetscBool sortIdx) 397 { 398 PetscInt v, w, vStart, vEnd, c, d; 399 PetscInt nVerts; 400 PetscInt closureSize = 0; 401 PetscInt *closure = NULL; 402 PetscInt *closureOrder; 403 PetscInt *invClosureOrder; 404 PetscInt *revlexOrder; 405 PetscInt *newNodeIdx; 406 PetscInt dim; 407 Vec coordVec; 408 const PetscScalar *coords; 409 PetscErrorCode ierr; 410 411 PetscFunctionBegin; 412 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 413 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 414 nVerts = vEnd - vStart; 415 ierr = PetscMalloc1(nVerts, &closureOrder);CHKERRQ(ierr); 416 ierr = PetscMalloc1(nVerts, &invClosureOrder);CHKERRQ(ierr); 417 ierr = PetscMalloc1(nVerts, &revlexOrder);CHKERRQ(ierr); 418 if (sortIdx) { /* bubble sort nodeIdx into revlex order */ 419 PetscInt nodeIdxDim = ni->nodeIdxDim; 420 PetscInt *idxOrder; 421 422 ierr = PetscMalloc1(nVerts * nodeIdxDim, &newNodeIdx);CHKERRQ(ierr); 423 ierr = PetscMalloc1(nVerts, &idxOrder);CHKERRQ(ierr); 424 for (v = 0; v < nVerts; v++) idxOrder[v] = v; 425 for (v = 0; v < nVerts; v++) { 426 for (w = v + 1; w < nVerts; w++) { 427 const PetscInt *iv = &(ni->nodeIdx[idxOrder[v] * nodeIdxDim]); 428 const PetscInt *iw = &(ni->nodeIdx[idxOrder[w] * nodeIdxDim]); 429 PetscInt diff = 0; 430 431 for (d = nodeIdxDim - 1; d >= 0; d--) if ((diff = (iv[d] - iw[d]))) break; 432 if (diff > 0) { 433 PetscInt swap = idxOrder[v]; 434 435 idxOrder[v] = idxOrder[w]; 436 idxOrder[w] = swap; 437 } 438 } 439 } 440 for (v = 0; v < nVerts; v++) { 441 for (d = 0; d < nodeIdxDim; d++) { 442 newNodeIdx[v * ni->nodeIdxDim + d] = ni->nodeIdx[idxOrder[v] * nodeIdxDim + d]; 443 } 444 } 445 ierr = PetscFree(ni->nodeIdx);CHKERRQ(ierr); 446 ni->nodeIdx = newNodeIdx; 447 newNodeIdx = NULL; 448 ierr = PetscFree(idxOrder);CHKERRQ(ierr); 449 } 450 ierr = DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 451 c = closureSize - nVerts; 452 for (v = 0; v < nVerts; v++) closureOrder[v] = closure[2 * (c + v)] - vStart; 453 for (v = 0; v < nVerts; v++) invClosureOrder[closureOrder[v]] = v; 454 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 455 ierr = DMGetCoordinatesLocal(dm, &coordVec);CHKERRQ(ierr); 456 ierr = VecGetArrayRead(coordVec, &coords);CHKERRQ(ierr); 457 /* bubble sort closure vertices by coordinates in revlex order */ 458 for (v = 0; v < nVerts; v++) revlexOrder[v] = v; 459 for (v = 0; v < nVerts; v++) { 460 for (w = v + 1; w < nVerts; w++) { 461 const PetscScalar *cv = &coords[closureOrder[revlexOrder[v]] * dim]; 462 const PetscScalar *cw = &coords[closureOrder[revlexOrder[w]] * dim]; 463 PetscReal diff = 0; 464 465 for (d = dim - 1; d >= 0; d--) if ((diff = PetscRealPart(cv[d] - cw[d])) != 0.) break; 466 if (diff > 0.) { 467 PetscInt swap = revlexOrder[v]; 468 469 revlexOrder[v] = revlexOrder[w]; 470 revlexOrder[w] = swap; 471 } 472 } 473 } 474 ierr = VecRestoreArrayRead(coordVec, &coords);CHKERRQ(ierr); 475 ierr = PetscMalloc1(ni->nodeIdxDim * nVerts, &newNodeIdx);CHKERRQ(ierr); 476 /* reorder nodeIdx to be in closure order */ 477 for (v = 0; v < nVerts; v++) { 478 for (d = 0; d < ni->nodeIdxDim; d++) { 479 newNodeIdx[revlexOrder[v] * ni->nodeIdxDim + d] = ni->nodeIdx[v * ni->nodeIdxDim + d]; 480 } 481 } 482 ierr = PetscFree(ni->nodeIdx);CHKERRQ(ierr); 483 ni->nodeIdx = newNodeIdx; 484 ni->perm = invClosureOrder; 485 ierr = PetscFree(revlexOrder);CHKERRQ(ierr); 486 ierr = PetscFree(closureOrder);CHKERRQ(ierr); 487 PetscFunctionReturn(0); 488 } 489 490 /* the coordinates of the simplex vertices are the corners of the barycentric simplex. 491 * When we stack them on top of each other in revlex order, they look like the identity matrix */ 492 static PetscErrorCode PetscLagNodeIndicesCreateSimplexVertices(DM dm, PetscLagNodeIndices *nodeIndices) 493 { 494 PetscLagNodeIndices ni; 495 PetscInt dim, d; 496 497 PetscErrorCode ierr; 498 499 PetscFunctionBegin; 500 ierr = PetscNew(&ni);CHKERRQ(ierr); 501 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 502 ni->nodeIdxDim = dim + 1; 503 ni->nodeVecDim = 0; 504 ni->nNodes = dim + 1; 505 ni->refct = 1; 506 ierr = PetscCalloc1((dim + 1)*(dim + 1), &(ni->nodeIdx));CHKERRQ(ierr); 507 for (d = 0; d < dim + 1; d++) ni->nodeIdx[d*(dim + 2)] = 1; 508 ierr = PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_FALSE);CHKERRQ(ierr); 509 *nodeIndices = ni; 510 PetscFunctionReturn(0); 511 } 512 513 /* A polytope that is a tensor product of a facet and a segment. 514 * We take whatever coordinate system was being used for the facet 515 * and we concatenate the barycentric coordinates for the vertices 516 * at the end of the segment, (1,0) and (0,1), to get a coordinate 517 * system for the tensor product element */ 518 static PetscErrorCode PetscLagNodeIndicesCreateTensorVertices(DM dm, PetscLagNodeIndices facetni, PetscLagNodeIndices *nodeIndices) 519 { 520 PetscLagNodeIndices ni; 521 PetscInt nodeIdxDim, subNodeIdxDim = facetni->nodeIdxDim; 522 PetscInt nVerts, nSubVerts = facetni->nNodes; 523 PetscInt dim, d, e, f, g; 524 525 PetscErrorCode ierr; 526 527 PetscFunctionBegin; 528 ierr = PetscNew(&ni);CHKERRQ(ierr); 529 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 530 ni->nodeIdxDim = nodeIdxDim = subNodeIdxDim + 2; 531 ni->nodeVecDim = 0; 532 ni->nNodes = nVerts = 2 * nSubVerts; 533 ni->refct = 1; 534 ierr = PetscCalloc1(nodeIdxDim * nVerts, &(ni->nodeIdx));CHKERRQ(ierr); 535 for (f = 0, d = 0; d < 2; d++) { 536 for (e = 0; e < nSubVerts; e++, f++) { 537 for (g = 0; g < subNodeIdxDim; g++) { 538 ni->nodeIdx[f * nodeIdxDim + g] = facetni->nodeIdx[e * subNodeIdxDim + g]; 539 } 540 ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim] = (1 - d); 541 ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim + 1] = d; 542 } 543 } 544 ierr = PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_TRUE);CHKERRQ(ierr); 545 *nodeIndices = ni; 546 PetscFunctionReturn(0); 547 } 548 549 /* This helps us compute symmetries, and it also helps us compute coordinates for dofs that are being pushed 550 * forward from a boundary mesh point. 551 * 552 * Input: 553 * 554 * dm - the target reference cell where we want new coordinates and dof directions to be valid 555 * vert - the vertex coordinate system for the target reference cell 556 * p - the point in the target reference cell that the dofs are coming from 557 * vertp - the vertex coordinate system for p's reference cell 558 * ornt - the resulting coordinates and dof vectors will be for p under this orientation 559 * nodep - the node coordinates and dof vectors in p's reference cell 560 * formDegree - the form degree that the dofs transform as 561 * 562 * Output: 563 * 564 * pfNodeIdx - the node coordinates for p's dofs, in the dm reference cell, from the ornt perspective 565 * pfNodeVec - the node dof vectors for p's dofs, in the dm reference cell, from the ornt perspective 566 */ 567 static PetscErrorCode PetscLagNodeIndicesPushForward(DM dm, PetscLagNodeIndices vert, PetscInt p, PetscLagNodeIndices vertp, PetscLagNodeIndices nodep, PetscInt ornt, PetscInt formDegree, PetscInt pfNodeIdx[], PetscReal pfNodeVec[]) 568 { 569 PetscInt *closureVerts; 570 PetscInt closureSize = 0; 571 PetscInt *closure = NULL; 572 PetscInt dim, pdim, c, i, j, k, n, v, vStart, vEnd; 573 PetscInt nSubVert = vertp->nNodes; 574 PetscInt nodeIdxDim = vert->nodeIdxDim; 575 PetscInt subNodeIdxDim = vertp->nodeIdxDim; 576 PetscInt nNodes = nodep->nNodes; 577 const PetscInt *vertIdx = vert->nodeIdx; 578 const PetscInt *subVertIdx = vertp->nodeIdx; 579 const PetscInt *nodeIdx = nodep->nodeIdx; 580 const PetscReal *nodeVec = nodep->nodeVec; 581 PetscReal *J, *Jstar; 582 PetscReal detJ; 583 PetscInt depth, pdepth, Nk, pNk; 584 Vec coordVec; 585 PetscScalar *newCoords = NULL; 586 const PetscScalar *oldCoords = NULL; 587 PetscErrorCode ierr; 588 589 PetscFunctionBegin; 590 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 591 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 592 ierr = DMGetCoordinatesLocal(dm, &coordVec);CHKERRQ(ierr); 593 ierr = DMPlexGetPointDepth(dm, p, &pdepth);CHKERRQ(ierr); 594 pdim = pdepth != depth ? pdepth != 0 ? pdepth : 0 : dim; 595 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 596 ierr = DMGetWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);CHKERRQ(ierr); 597 ierr = DMPlexGetTransitiveClosure_Internal(dm, p, ornt, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 598 c = closureSize - nSubVert; 599 /* we want which cell closure indices the closure of this point corresponds to */ 600 for (v = 0; v < nSubVert; v++) closureVerts[v] = vert->perm[closure[2 * (c + v)] - vStart]; 601 ierr = DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 602 /* push forward indices */ 603 for (i = 0; i < nodeIdxDim; i++) { /* for every component of the target index space */ 604 /* check if this is a component that all vertices around this point have in common */ 605 for (j = 1; j < nSubVert; j++) { 606 if (vertIdx[closureVerts[j] * nodeIdxDim + i] != vertIdx[closureVerts[0] * nodeIdxDim + i]) break; 607 } 608 if (j == nSubVert) { /* all vertices have this component in common, directly copy to output */ 609 PetscInt val = vertIdx[closureVerts[0] * nodeIdxDim + i]; 610 for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = val; 611 } else { 612 PetscInt subi = -1; 613 /* there must be a component in vertp that looks the same */ 614 for (k = 0; k < subNodeIdxDim; k++) { 615 for (j = 0; j < nSubVert; j++) { 616 if (vertIdx[closureVerts[j] * nodeIdxDim + i] != subVertIdx[j * subNodeIdxDim + k]) break; 617 } 618 if (j == nSubVert) { 619 subi = k; 620 break; 621 } 622 } 623 if (subi < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not find matching coordinate\n"); 624 /* that component in the vertp system becomes component i in the vert system for each dof */ 625 for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = nodeIdx[n * subNodeIdxDim + subi]; 626 } 627 } 628 /* push forward vectors */ 629 ierr = DMGetWorkArray(dm, dim * dim, MPIU_REAL, &J);CHKERRQ(ierr); 630 if (ornt != 0) { /* temporarily change the coordinate vector so 631 DMPlexComputeCellGeometryAffineFEM gives us the Jacobian we want */ 632 PetscInt closureSize2 = 0; 633 PetscInt *closure2 = NULL; 634 635 ierr = DMPlexGetTransitiveClosure_Internal(dm, p, 0, PETSC_TRUE, &closureSize2, &closure2);CHKERRQ(ierr); 636 ierr = PetscMalloc1(dim * nSubVert, &newCoords);CHKERRQ(ierr); 637 ierr = VecGetArrayRead(coordVec, &oldCoords);CHKERRQ(ierr); 638 for (v = 0; v < nSubVert; v++) { 639 PetscInt d; 640 for (d = 0; d < dim; d++) { 641 newCoords[(closure2[2 * (c + v)] - vStart) * dim + d] = oldCoords[closureVerts[v] * dim + d]; 642 } 643 } 644 ierr = VecRestoreArrayRead(coordVec, &oldCoords);CHKERRQ(ierr); 645 ierr = DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize2, &closure2);CHKERRQ(ierr); 646 ierr = VecPlaceArray(coordVec, newCoords);CHKERRQ(ierr); 647 } 648 ierr = DMPlexComputeCellGeometryAffineFEM(dm, p, NULL, J, NULL, &detJ);CHKERRQ(ierr); 649 if (ornt != 0) { 650 ierr = VecResetArray(coordVec);CHKERRQ(ierr); 651 ierr = PetscFree(newCoords);CHKERRQ(ierr); 652 } 653 ierr = DMRestoreWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);CHKERRQ(ierr); 654 /* compactify */ 655 for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j]; 656 /* We have the Jacobian mapping the point's reference cell to this reference cell: 657 * pulling back a function to the point and applying the dof is what we want, 658 * so we get the pullback matrix and multiply the dof by that matrix on the right */ 659 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);CHKERRQ(ierr); 660 ierr = PetscDTBinomialInt(pdim, PetscAbsInt(formDegree), &pNk);CHKERRQ(ierr); 661 ierr = DMGetWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);CHKERRQ(ierr); 662 ierr = PetscDTAltVPullbackMatrix(pdim, dim, J, formDegree, Jstar);CHKERRQ(ierr); 663 for (n = 0; n < nNodes; n++) { 664 for (i = 0; i < Nk; i++) { 665 PetscReal val = 0.; 666 for (j = 0; j < pNk; j++) val += nodeVec[n * pNk + j] * Jstar[j * Nk + i]; 667 pfNodeVec[n * Nk + i] = val; 668 } 669 } 670 ierr = DMRestoreWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);CHKERRQ(ierr); 671 ierr = DMRestoreWorkArray(dm, dim * dim, MPIU_REAL, &J);CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 /* given to sets of nodes, take the tensor product, where the product of the dof indices is concatenation and the 676 * product of the dof vectors is the wedge product */ 677 static PetscErrorCode PetscLagNodeIndicesTensor(PetscLagNodeIndices tracei, PetscInt dimT, PetscInt kT, PetscLagNodeIndices fiberi, PetscInt dimF, PetscInt kF, PetscLagNodeIndices *nodeIndices) 678 { 679 PetscInt dim = dimT + dimF; 680 PetscInt nodeIdxDim, nNodes; 681 PetscInt formDegree = kT + kF; 682 PetscInt Nk, NkT, NkF; 683 PetscInt MkT, MkF; 684 PetscLagNodeIndices ni; 685 PetscInt i, j, l; 686 PetscReal *projF, *projT; 687 PetscReal *projFstar, *projTstar; 688 PetscReal *workF, *workF2, *workT, *workT2, *work, *work2; 689 PetscReal *wedgeMat; 690 PetscReal sign; 691 PetscErrorCode ierr; 692 693 PetscFunctionBegin; 694 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);CHKERRQ(ierr); 695 ierr = PetscDTBinomialInt(dimT, PetscAbsInt(kT), &NkT);CHKERRQ(ierr); 696 ierr = PetscDTBinomialInt(dimF, PetscAbsInt(kF), &NkF);CHKERRQ(ierr); 697 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kT), &MkT);CHKERRQ(ierr); 698 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kF), &MkF);CHKERRQ(ierr); 699 ierr = PetscNew(&ni);CHKERRQ(ierr); 700 ni->nodeIdxDim = nodeIdxDim = tracei->nodeIdxDim + fiberi->nodeIdxDim; 701 ni->nodeVecDim = Nk; 702 ni->nNodes = nNodes = tracei->nNodes * fiberi->nNodes; 703 ni->refct = 1; 704 ierr = PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));CHKERRQ(ierr); 705 /* first concatenate the indices */ 706 for (l = 0, j = 0; j < fiberi->nNodes; j++) { 707 for (i = 0; i < tracei->nNodes; i++, l++) { 708 PetscInt m, n = 0; 709 710 for (m = 0; m < tracei->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = tracei->nodeIdx[i * tracei->nodeIdxDim + m]; 711 for (m = 0; m < fiberi->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = fiberi->nodeIdx[j * fiberi->nodeIdxDim + m]; 712 } 713 } 714 715 /* now wedge together the push-forward vectors */ 716 ierr = PetscMalloc1(nNodes * Nk, &(ni->nodeVec));CHKERRQ(ierr); 717 ierr = PetscCalloc2(dimT*dim, &projT, dimF*dim, &projF);CHKERRQ(ierr); 718 for (i = 0; i < dimT; i++) projT[i * (dim + 1)] = 1.; 719 for (i = 0; i < dimF; i++) projF[i * (dim + dimT + 1) + dimT] = 1.; 720 ierr = PetscMalloc2(MkT*NkT, &projTstar, MkF*NkF, &projFstar);CHKERRQ(ierr); 721 ierr = PetscDTAltVPullbackMatrix(dim, dimT, projT, kT, projTstar);CHKERRQ(ierr); 722 ierr = PetscDTAltVPullbackMatrix(dim, dimF, projF, kF, projFstar);CHKERRQ(ierr); 723 ierr = PetscMalloc6(MkT, &workT, MkT, &workT2, MkF, &workF, MkF, &workF2, Nk, &work, Nk, &work2);CHKERRQ(ierr); 724 ierr = PetscMalloc1(Nk * MkT, &wedgeMat);CHKERRQ(ierr); 725 sign = (PetscAbsInt(kT * kF) & 1) ? -1. : 1.; 726 for (l = 0, j = 0; j < fiberi->nNodes; j++) { 727 PetscInt d, e; 728 729 /* push forward fiber k-form */ 730 for (d = 0; d < MkF; d++) { 731 PetscReal val = 0.; 732 for (e = 0; e < NkF; e++) val += projFstar[d * NkF + e] * fiberi->nodeVec[j * NkF + e]; 733 workF[d] = val; 734 } 735 /* Hodge star to proper form if necessary */ 736 if (kF < 0) { 737 for (d = 0; d < MkF; d++) workF2[d] = workF[d]; 738 ierr = PetscDTAltVStar(dim, PetscAbsInt(kF), 1, workF2, workF);CHKERRQ(ierr); 739 } 740 /* Compute the matrix that wedges this form with one of the trace k-form */ 741 ierr = PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kF), PetscAbsInt(kT), workF, wedgeMat);CHKERRQ(ierr); 742 for (i = 0; i < tracei->nNodes; i++, l++) { 743 /* push forward trace k-form */ 744 for (d = 0; d < MkT; d++) { 745 PetscReal val = 0.; 746 for (e = 0; e < NkT; e++) val += projTstar[d * NkT + e] * tracei->nodeVec[i * NkT + e]; 747 workT[d] = val; 748 } 749 /* Hodge star to proper form if necessary */ 750 if (kT < 0) { 751 for (d = 0; d < MkT; d++) workT2[d] = workT[d]; 752 ierr = PetscDTAltVStar(dim, PetscAbsInt(kT), 1, workT2, workT);CHKERRQ(ierr); 753 } 754 /* compute the wedge product of the push-forward trace form and firer forms */ 755 for (d = 0; d < Nk; d++) { 756 PetscReal val = 0.; 757 for (e = 0; e < MkT; e++) val += wedgeMat[d * MkT + e] * workT[e]; 758 work[d] = val; 759 } 760 /* inverse Hodge star from proper form if necessary */ 761 if (formDegree < 0) { 762 for (d = 0; d < Nk; d++) work2[d] = work[d]; 763 ierr = PetscDTAltVStar(dim, PetscAbsInt(formDegree), -1, work2, work);CHKERRQ(ierr); 764 } 765 /* insert into the array (adjusting for sign) */ 766 for (d = 0; d < Nk; d++) ni->nodeVec[l * Nk + d] = sign * work[d]; 767 } 768 } 769 ierr = PetscFree(wedgeMat);CHKERRQ(ierr); 770 ierr = PetscFree6(workT, workT2, workF, workF2, work, work2);CHKERRQ(ierr); 771 ierr = PetscFree2(projTstar, projFstar);CHKERRQ(ierr); 772 ierr = PetscFree2(projT, projF);CHKERRQ(ierr); 773 *nodeIndices = ni; 774 PetscFunctionReturn(0); 775 } 776 777 /* simple union of two sets of nodes */ 778 static PetscErrorCode PetscLagNodeIndicesMerge(PetscLagNodeIndices niA, PetscLagNodeIndices niB, PetscLagNodeIndices *nodeIndices) 779 { 780 PetscLagNodeIndices ni; 781 PetscInt nodeIdxDim, nodeVecDim, nNodes; 782 PetscErrorCode ierr; 783 784 PetscFunctionBegin; 785 ierr = PetscNew(&ni);CHKERRQ(ierr); 786 ni->nodeIdxDim = nodeIdxDim = niA->nodeIdxDim; 787 if (niB->nodeIdxDim != nodeIdxDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeIdxDim"); 788 ni->nodeVecDim = nodeVecDim = niA->nodeVecDim; 789 if (niB->nodeVecDim != nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeVecDim"); 790 ni->nNodes = nNodes = niA->nNodes + niB->nNodes; 791 ni->refct = 1; 792 ierr = PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));CHKERRQ(ierr); 793 ierr = PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));CHKERRQ(ierr); 794 ierr = PetscArraycpy(ni->nodeIdx, niA->nodeIdx, niA->nNodes * nodeIdxDim);CHKERRQ(ierr); 795 ierr = PetscArraycpy(ni->nodeVec, niA->nodeVec, niA->nNodes * nodeVecDim);CHKERRQ(ierr); 796 ierr = PetscArraycpy(&(ni->nodeIdx[niA->nNodes * nodeIdxDim]), niB->nodeIdx, niB->nNodes * nodeIdxDim);CHKERRQ(ierr); 797 ierr = PetscArraycpy(&(ni->nodeVec[niA->nNodes * nodeVecDim]), niB->nodeVec, niB->nNodes * nodeVecDim);CHKERRQ(ierr); 798 *nodeIndices = ni; 799 PetscFunctionReturn(0); 800 } 801 802 #define PETSCTUPINTCOMPREVLEX(N) \ 803 static int PetscTupIntCompRevlex_##N(const void *a, const void *b) \ 804 { \ 805 const PetscInt *A = (const PetscInt *) a; \ 806 const PetscInt *B = (const PetscInt *) b; \ 807 int i; \ 808 PetscInt diff = 0; \ 809 for (i = 0; i < N; i++) { \ 810 diff = A[N - i] - B[N - i]; \ 811 if (diff) break; \ 812 } \ 813 return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1; \ 814 } 815 816 PETSCTUPINTCOMPREVLEX(3) 817 PETSCTUPINTCOMPREVLEX(4) 818 PETSCTUPINTCOMPREVLEX(5) 819 PETSCTUPINTCOMPREVLEX(6) 820 PETSCTUPINTCOMPREVLEX(7) 821 822 static int PetscTupIntCompRevlex_N(const void *a, const void *b) 823 { 824 const PetscInt *A = (const PetscInt *) a; 825 const PetscInt *B = (const PetscInt *) b; 826 int i; 827 int N = A[0]; 828 PetscInt diff = 0; 829 for (i = 0; i < N; i++) { 830 diff = A[N - i] - B[N - i]; 831 if (diff) break; 832 } 833 return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1; 834 } 835 836 /* The nodes are not necessarily in revlex order wrt nodeIdx: get the permutation 837 * that puts them in that order */ 838 static PetscErrorCode PetscLagNodeIndicesGetPermutation(PetscLagNodeIndices ni, PetscInt *perm[]) 839 { 840 PetscErrorCode ierr; 841 842 PetscFunctionBegin; 843 if (!(ni->perm)) { 844 PetscInt *sorter; 845 PetscInt m = ni->nNodes; 846 PetscInt nodeIdxDim = ni->nodeIdxDim; 847 PetscInt i, j, k, l; 848 PetscInt *prm; 849 int (*comp) (const void *, const void *); 850 851 ierr = PetscMalloc1((nodeIdxDim + 2) * m, &sorter);CHKERRQ(ierr); 852 for (k = 0, l = 0, i = 0; i < m; i++) { 853 sorter[k++] = nodeIdxDim + 1; 854 sorter[k++] = i; 855 for (j = 0; j < nodeIdxDim; j++) sorter[k++] = ni->nodeIdx[l++]; 856 } 857 switch (nodeIdxDim) { 858 case 2: 859 comp = PetscTupIntCompRevlex_3; 860 break; 861 case 3: 862 comp = PetscTupIntCompRevlex_4; 863 break; 864 case 4: 865 comp = PetscTupIntCompRevlex_5; 866 break; 867 case 5: 868 comp = PetscTupIntCompRevlex_6; 869 break; 870 case 6: 871 comp = PetscTupIntCompRevlex_7; 872 break; 873 default: 874 comp = PetscTupIntCompRevlex_N; 875 break; 876 } 877 qsort(sorter, m, (nodeIdxDim + 2) * sizeof(PetscInt), comp); 878 ierr = PetscMalloc1(m, &prm);CHKERRQ(ierr); 879 for (i = 0; i < m; i++) prm[i] = sorter[(nodeIdxDim + 2) * i + 1]; 880 ni->perm = prm; 881 ierr = PetscFree(sorter);CHKERRQ(ierr); 882 } 883 *perm = ni->perm; 884 PetscFunctionReturn(0); 885 } 886 887 static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp) 888 { 889 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 890 PetscErrorCode ierr; 891 892 PetscFunctionBegin; 893 if (lag->symperms) { 894 PetscInt **selfSyms = lag->symperms[0]; 895 896 if (selfSyms) { 897 PetscInt i, **allocated = &selfSyms[-lag->selfSymOff]; 898 899 for (i = 0; i < lag->numSelfSym; i++) { 900 ierr = PetscFree(allocated[i]);CHKERRQ(ierr); 901 } 902 ierr = PetscFree(allocated);CHKERRQ(ierr); 903 } 904 ierr = PetscFree(lag->symperms);CHKERRQ(ierr); 905 } 906 if (lag->symflips) { 907 PetscScalar **selfSyms = lag->symflips[0]; 908 909 if (selfSyms) { 910 PetscInt i; 911 PetscScalar **allocated = &selfSyms[-lag->selfSymOff]; 912 913 for (i = 0; i < lag->numSelfSym; i++) { 914 ierr = PetscFree(allocated[i]);CHKERRQ(ierr); 915 } 916 ierr = PetscFree(allocated);CHKERRQ(ierr); 917 } 918 ierr = PetscFree(lag->symflips);CHKERRQ(ierr); 919 } 920 ierr = Petsc1DNodeFamilyDestroy(&(lag->nodeFamily));CHKERRQ(ierr); 921 ierr = PetscLagNodeIndicesDestroy(&(lag->vertIndices));CHKERRQ(ierr); 922 ierr = PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));CHKERRQ(ierr); 923 ierr = PetscLagNodeIndicesDestroy(&(lag->allNodeIndices));CHKERRQ(ierr); 924 ierr = PetscFree(lag);CHKERRQ(ierr); 925 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);CHKERRQ(ierr); 926 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);CHKERRQ(ierr); 927 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);CHKERRQ(ierr); 928 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);CHKERRQ(ierr); 929 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL);CHKERRQ(ierr); 930 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL);CHKERRQ(ierr); 931 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL);CHKERRQ(ierr); 932 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL);CHKERRQ(ierr); 933 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetUseMoments_C", NULL);CHKERRQ(ierr); 934 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetUseMoments_C", NULL);CHKERRQ(ierr); 935 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetMomentOrder_C", NULL);CHKERRQ(ierr); 936 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetMomentOrder_C", NULL);CHKERRQ(ierr); 937 PetscFunctionReturn(0); 938 } 939 940 static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer) 941 { 942 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 943 PetscErrorCode ierr; 944 945 PetscFunctionBegin; 946 ierr = PetscViewerASCIIPrintf(viewer, "%s %s%sLagrange dual space\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "", lag->trimmed ? "trimmed " : "");CHKERRQ(ierr); 947 PetscFunctionReturn(0); 948 } 949 950 static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer) 951 { 952 PetscBool iascii; 953 PetscErrorCode ierr; 954 955 PetscFunctionBegin; 956 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 957 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 958 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 959 if (iascii) {ierr = PetscDualSpaceLagrangeView_Ascii(sp, viewer);CHKERRQ(ierr);} 960 PetscFunctionReturn(0); 961 } 962 963 static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp) 964 { 965 PetscBool continuous, tensor, trimmed, flg, flg2, flg3; 966 PetscDTNodeType nodeType; 967 PetscReal nodeExponent; 968 PetscInt momentOrder; 969 PetscBool nodeEndpoints, useMoments; 970 PetscErrorCode ierr; 971 972 PetscFunctionBegin; 973 ierr = PetscDualSpaceLagrangeGetContinuity(sp, &continuous);CHKERRQ(ierr); 974 ierr = PetscDualSpaceLagrangeGetTensor(sp, &tensor);CHKERRQ(ierr); 975 ierr = PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);CHKERRQ(ierr); 976 ierr = PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &nodeEndpoints, &nodeExponent);CHKERRQ(ierr); 977 if (nodeType == PETSCDTNODES_DEFAULT) nodeType = PETSCDTNODES_GAUSSJACOBI; 978 ierr = PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments);CHKERRQ(ierr); 979 ierr = PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder);CHKERRQ(ierr); 980 ierr = PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");CHKERRQ(ierr); 981 ierr = PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);CHKERRQ(ierr); 982 if (flg) {ierr = PetscDualSpaceLagrangeSetContinuity(sp, continuous);CHKERRQ(ierr);} 983 ierr = PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetTensor", tensor, &tensor, &flg);CHKERRQ(ierr); 984 if (flg) {ierr = PetscDualSpaceLagrangeSetTensor(sp, tensor);CHKERRQ(ierr);} 985 ierr = PetscOptionsBool("-petscdualspace_lagrange_trimmed", "Flag for trimmed dual space", "PetscDualSpaceLagrangeSetTrimmed", trimmed, &trimmed, &flg);CHKERRQ(ierr); 986 if (flg) {ierr = PetscDualSpaceLagrangeSetTrimmed(sp, trimmed);CHKERRQ(ierr);} 987 ierr = PetscOptionsEnum("-petscdualspace_lagrange_node_type", "Lagrange node location type", "PetscDualSpaceLagrangeSetNodeType", PetscDTNodeTypes, (PetscEnum)nodeType, (PetscEnum *)&nodeType, &flg);CHKERRQ(ierr); 988 ierr = PetscOptionsBool("-petscdualspace_lagrange_node_endpoints", "Flag for nodes that include endpoints", "PetscDualSpaceLagrangeSetNodeType", nodeEndpoints, &nodeEndpoints, &flg2);CHKERRQ(ierr); 989 flg3 = PETSC_FALSE; 990 if (nodeType == PETSCDTNODES_GAUSSJACOBI) { 991 ierr = PetscOptionsReal("-petscdualspace_lagrange_node_exponent", "Gauss-Jacobi weight function exponent", "PetscDualSpaceLagrangeSetNodeType", nodeExponent, &nodeExponent, &flg3);CHKERRQ(ierr); 992 } 993 if (flg || flg2 || flg3) {ierr = PetscDualSpaceLagrangeSetNodeType(sp, nodeType, nodeEndpoints, nodeExponent);CHKERRQ(ierr);} 994 ierr = PetscOptionsBool("-petscdualspace_lagrange_use_moments", "Use moments (where appropriate) for functionals", "PetscDualSpaceLagrangeSetUseMoments", useMoments, &useMoments, &flg);CHKERRQ(ierr); 995 if (flg) {ierr = PetscDualSpaceLagrangeSetUseMoments(sp, useMoments);CHKERRQ(ierr);} 996 ierr = PetscOptionsInt("-petscdualspace_lagrange_moment_order", "Quadrature order for moment functionals", "PetscDualSpaceLagrangeSetMomentOrder", momentOrder, &momentOrder, &flg);CHKERRQ(ierr); 997 if (flg) {ierr = PetscDualSpaceLagrangeSetMomentOrder(sp, momentOrder);CHKERRQ(ierr);} 998 ierr = PetscOptionsTail();CHKERRQ(ierr); 999 PetscFunctionReturn(0); 1000 } 1001 1002 static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace spNew) 1003 { 1004 PetscBool cont, tensor, trimmed, boundary; 1005 PetscDTNodeType nodeType; 1006 PetscReal exponent; 1007 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 1008 PetscErrorCode ierr; 1009 1010 PetscFunctionBegin; 1011 ierr = PetscDualSpaceLagrangeGetContinuity(sp, &cont);CHKERRQ(ierr); 1012 ierr = PetscDualSpaceLagrangeSetContinuity(spNew, cont);CHKERRQ(ierr); 1013 ierr = PetscDualSpaceLagrangeGetTensor(sp, &tensor);CHKERRQ(ierr); 1014 ierr = PetscDualSpaceLagrangeSetTensor(spNew, tensor);CHKERRQ(ierr); 1015 ierr = PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);CHKERRQ(ierr); 1016 ierr = PetscDualSpaceLagrangeSetTrimmed(spNew, trimmed);CHKERRQ(ierr); 1017 ierr = PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &boundary, &exponent);CHKERRQ(ierr); 1018 ierr = PetscDualSpaceLagrangeSetNodeType(spNew, nodeType, boundary, exponent);CHKERRQ(ierr); 1019 if (lag->nodeFamily) { 1020 PetscDualSpace_Lag *lagnew = (PetscDualSpace_Lag *) spNew->data; 1021 1022 ierr = Petsc1DNodeFamilyReference(lag->nodeFamily);CHKERRQ(ierr); 1023 lagnew->nodeFamily = lag->nodeFamily; 1024 } 1025 PetscFunctionReturn(0); 1026 } 1027 1028 /* for making tensor product spaces: take a dual space and product a segment space that has all the same 1029 * specifications (trimmed, continuous, order, node set), except for the form degree */ 1030 static PetscErrorCode PetscDualSpaceCreateEdgeSubspace_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt k, PetscInt Nc, PetscBool interiorOnly, PetscDualSpace *bdsp) 1031 { 1032 DM K; 1033 PetscDualSpace_Lag *newlag; 1034 PetscErrorCode ierr; 1035 1036 PetscFunctionBegin; 1037 ierr = PetscDualSpaceDuplicate(sp,bdsp);CHKERRQ(ierr); 1038 ierr = PetscDualSpaceSetFormDegree(*bdsp, k);CHKERRQ(ierr); 1039 ierr = PetscDualSpaceCreateReferenceCell(*bdsp, 1, PETSC_TRUE, &K);CHKERRQ(ierr); 1040 ierr = PetscDualSpaceSetDM(*bdsp, K);CHKERRQ(ierr); 1041 ierr = DMDestroy(&K);CHKERRQ(ierr); 1042 ierr = PetscDualSpaceSetOrder(*bdsp, order);CHKERRQ(ierr); 1043 ierr = PetscDualSpaceSetNumComponents(*bdsp, Nc);CHKERRQ(ierr); 1044 newlag = (PetscDualSpace_Lag *) (*bdsp)->data; 1045 newlag->interiorOnly = interiorOnly; 1046 ierr = PetscDualSpaceSetUp(*bdsp);CHKERRQ(ierr); 1047 PetscFunctionReturn(0); 1048 } 1049 1050 /* just the points, weights aren't handled */ 1051 static PetscErrorCode PetscQuadratureCreateTensor(PetscQuadrature trace, PetscQuadrature fiber, PetscQuadrature *product) 1052 { 1053 PetscInt dimTrace, dimFiber; 1054 PetscInt numPointsTrace, numPointsFiber; 1055 PetscInt dim, numPoints; 1056 const PetscReal *pointsTrace; 1057 const PetscReal *pointsFiber; 1058 PetscReal *points; 1059 PetscInt i, j, k, p; 1060 PetscErrorCode ierr; 1061 1062 PetscFunctionBegin; 1063 ierr = PetscQuadratureGetData(trace, &dimTrace, NULL, &numPointsTrace, &pointsTrace, NULL);CHKERRQ(ierr); 1064 ierr = PetscQuadratureGetData(fiber, &dimFiber, NULL, &numPointsFiber, &pointsFiber, NULL);CHKERRQ(ierr); 1065 dim = dimTrace + dimFiber; 1066 numPoints = numPointsFiber * numPointsTrace; 1067 ierr = PetscMalloc1(numPoints * dim, &points);CHKERRQ(ierr); 1068 for (p = 0, j = 0; j < numPointsFiber; j++) { 1069 for (i = 0; i < numPointsTrace; i++, p++) { 1070 for (k = 0; k < dimTrace; k++) points[p * dim + k] = pointsTrace[i * dimTrace + k]; 1071 for (k = 0; k < dimFiber; k++) points[p * dim + dimTrace + k] = pointsFiber[j * dimFiber + k]; 1072 } 1073 } 1074 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, product);CHKERRQ(ierr); 1075 ierr = PetscQuadratureSetData(*product, dim, 0, numPoints, points, NULL);CHKERRQ(ierr); 1076 PetscFunctionReturn(0); 1077 } 1078 1079 /* Kronecker tensor product where matrix is considered a matrix of k-forms, so that 1080 * the entries in the product matrix are wedge products of the entries in the original matrices */ 1081 static PetscErrorCode MatTensorAltV(Mat trace, Mat fiber, PetscInt dimTrace, PetscInt kTrace, PetscInt dimFiber, PetscInt kFiber, Mat *product) 1082 { 1083 PetscInt mTrace, nTrace, mFiber, nFiber, m, n, k, i, j, l; 1084 PetscInt dim, NkTrace, NkFiber, Nk; 1085 PetscInt dT, dF; 1086 PetscInt *nnzTrace, *nnzFiber, *nnz; 1087 PetscInt iT, iF, jT, jF, il, jl; 1088 PetscReal *workT, *workT2, *workF, *workF2, *work, *workstar; 1089 PetscReal *projT, *projF; 1090 PetscReal *projTstar, *projFstar; 1091 PetscReal *wedgeMat; 1092 PetscReal sign; 1093 PetscScalar *workS; 1094 Mat prod; 1095 /* this produces dof groups that look like the identity */ 1096 PetscErrorCode ierr; 1097 1098 PetscFunctionBegin; 1099 ierr = MatGetSize(trace, &mTrace, &nTrace);CHKERRQ(ierr); 1100 ierr = PetscDTBinomialInt(dimTrace, PetscAbsInt(kTrace), &NkTrace);CHKERRQ(ierr); 1101 if (nTrace % NkTrace) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point value space of trace matrix is not a multiple of k-form size"); 1102 ierr = MatGetSize(fiber, &mFiber, &nFiber);CHKERRQ(ierr); 1103 ierr = PetscDTBinomialInt(dimFiber, PetscAbsInt(kFiber), &NkFiber);CHKERRQ(ierr); 1104 if (nFiber % NkFiber) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point value space of fiber matrix is not a multiple of k-form size"); 1105 ierr = PetscMalloc2(mTrace, &nnzTrace, mFiber, &nnzFiber);CHKERRQ(ierr); 1106 for (i = 0; i < mTrace; i++) { 1107 ierr = MatGetRow(trace, i, &(nnzTrace[i]), NULL, NULL);CHKERRQ(ierr); 1108 if (nnzTrace[i] % NkTrace) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in trace matrix are not in k-form size blocks"); 1109 } 1110 for (i = 0; i < mFiber; i++) { 1111 ierr = MatGetRow(fiber, i, &(nnzFiber[i]), NULL, NULL);CHKERRQ(ierr); 1112 if (nnzFiber[i] % NkFiber) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in fiber matrix are not in k-form size blocks"); 1113 } 1114 dim = dimTrace + dimFiber; 1115 k = kFiber + kTrace; 1116 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 1117 m = mTrace * mFiber; 1118 ierr = PetscMalloc1(m, &nnz);CHKERRQ(ierr); 1119 for (l = 0, j = 0; j < mFiber; j++) for (i = 0; i < mTrace; i++, l++) nnz[l] = (nnzTrace[i] / NkTrace) * (nnzFiber[j] / NkFiber) * Nk; 1120 n = (nTrace / NkTrace) * (nFiber / NkFiber) * Nk; 1121 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &prod);CHKERRQ(ierr); 1122 ierr = PetscFree(nnz);CHKERRQ(ierr); 1123 ierr = PetscFree2(nnzTrace,nnzFiber);CHKERRQ(ierr); 1124 /* reasoning about which points each dof needs depends on having zeros computed at points preserved */ 1125 ierr = MatSetOption(prod, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr); 1126 /* compute pullbacks */ 1127 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kTrace), &dT);CHKERRQ(ierr); 1128 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kFiber), &dF);CHKERRQ(ierr); 1129 ierr = PetscMalloc4(dimTrace * dim, &projT, dimFiber * dim, &projF, dT * NkTrace, &projTstar, dF * NkFiber, &projFstar);CHKERRQ(ierr); 1130 ierr = PetscArrayzero(projT, dimTrace * dim);CHKERRQ(ierr); 1131 for (i = 0; i < dimTrace; i++) projT[i * (dim + 1)] = 1.; 1132 ierr = PetscArrayzero(projF, dimFiber * dim);CHKERRQ(ierr); 1133 for (i = 0; i < dimFiber; i++) projF[i * (dim + 1) + dimTrace] = 1.; 1134 ierr = PetscDTAltVPullbackMatrix(dim, dimTrace, projT, kTrace, projTstar);CHKERRQ(ierr); 1135 ierr = PetscDTAltVPullbackMatrix(dim, dimFiber, projF, kFiber, projFstar);CHKERRQ(ierr); 1136 ierr = PetscMalloc5(dT, &workT, dF, &workF, Nk, &work, Nk, &workstar, Nk, &workS);CHKERRQ(ierr); 1137 ierr = PetscMalloc2(dT, &workT2, dF, &workF2);CHKERRQ(ierr); 1138 ierr = PetscMalloc1(Nk * dT, &wedgeMat);CHKERRQ(ierr); 1139 sign = (PetscAbsInt(kTrace * kFiber) & 1) ? -1. : 1.; 1140 for (i = 0, iF = 0; iF < mFiber; iF++) { 1141 PetscInt ncolsF, nformsF; 1142 const PetscInt *colsF; 1143 const PetscScalar *valsF; 1144 1145 ierr = MatGetRow(fiber, iF, &ncolsF, &colsF, &valsF);CHKERRQ(ierr); 1146 nformsF = ncolsF / NkFiber; 1147 for (iT = 0; iT < mTrace; iT++, i++) { 1148 PetscInt ncolsT, nformsT; 1149 const PetscInt *colsT; 1150 const PetscScalar *valsT; 1151 1152 ierr = MatGetRow(trace, iT, &ncolsT, &colsT, &valsT);CHKERRQ(ierr); 1153 nformsT = ncolsT / NkTrace; 1154 for (j = 0, jF = 0; jF < nformsF; jF++) { 1155 PetscInt colF = colsF[jF * NkFiber] / NkFiber; 1156 1157 for (il = 0; il < dF; il++) { 1158 PetscReal val = 0.; 1159 for (jl = 0; jl < NkFiber; jl++) val += projFstar[il * NkFiber + jl] * PetscRealPart(valsF[jF * NkFiber + jl]); 1160 workF[il] = val; 1161 } 1162 if (kFiber < 0) { 1163 for (il = 0; il < dF; il++) workF2[il] = workF[il]; 1164 ierr = PetscDTAltVStar(dim, PetscAbsInt(kFiber), 1, workF2, workF);CHKERRQ(ierr); 1165 } 1166 ierr = PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kFiber), PetscAbsInt(kTrace), workF, wedgeMat);CHKERRQ(ierr); 1167 for (jT = 0; jT < nformsT; jT++, j++) { 1168 PetscInt colT = colsT[jT * NkTrace] / NkTrace; 1169 PetscInt col = colF * (nTrace / NkTrace) + colT; 1170 const PetscScalar *vals; 1171 1172 for (il = 0; il < dT; il++) { 1173 PetscReal val = 0.; 1174 for (jl = 0; jl < NkTrace; jl++) val += projTstar[il * NkTrace + jl] * PetscRealPart(valsT[jT * NkTrace + jl]); 1175 workT[il] = val; 1176 } 1177 if (kTrace < 0) { 1178 for (il = 0; il < dT; il++) workT2[il] = workT[il]; 1179 ierr = PetscDTAltVStar(dim, PetscAbsInt(kTrace), 1, workT2, workT);CHKERRQ(ierr); 1180 } 1181 1182 for (il = 0; il < Nk; il++) { 1183 PetscReal val = 0.; 1184 for (jl = 0; jl < dT; jl++) val += sign * wedgeMat[il * dT + jl] * workT[jl]; 1185 work[il] = val; 1186 } 1187 if (k < 0) { 1188 ierr = PetscDTAltVStar(dim, PetscAbsInt(k), -1, work, workstar);CHKERRQ(ierr); 1189 #if defined(PETSC_USE_COMPLEX) 1190 for (l = 0; l < Nk; l++) workS[l] = workstar[l]; 1191 vals = &workS[0]; 1192 #else 1193 vals = &workstar[0]; 1194 #endif 1195 } else { 1196 #if defined(PETSC_USE_COMPLEX) 1197 for (l = 0; l < Nk; l++) workS[l] = work[l]; 1198 vals = &workS[0]; 1199 #else 1200 vals = &work[0]; 1201 #endif 1202 } 1203 for (l = 0; l < Nk; l++) { 1204 ierr = MatSetValue(prod, i, col * Nk + l, vals[l], INSERT_VALUES);CHKERRQ(ierr); 1205 } /* Nk */ 1206 } /* jT */ 1207 } /* jF */ 1208 ierr = MatRestoreRow(trace, iT, &ncolsT, &colsT, &valsT);CHKERRQ(ierr); 1209 } /* iT */ 1210 ierr = MatRestoreRow(fiber, iF, &ncolsF, &colsF, &valsF);CHKERRQ(ierr); 1211 } /* iF */ 1212 ierr = PetscFree(wedgeMat);CHKERRQ(ierr); 1213 ierr = PetscFree4(projT, projF, projTstar, projFstar);CHKERRQ(ierr); 1214 ierr = PetscFree2(workT2, workF2);CHKERRQ(ierr); 1215 ierr = PetscFree5(workT, workF, work, workstar, workS);CHKERRQ(ierr); 1216 ierr = MatAssemblyBegin(prod, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1217 ierr = MatAssemblyEnd(prod, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1218 *product = prod; 1219 PetscFunctionReturn(0); 1220 } 1221 1222 /* Union of quadrature points, with an attempt to identify commont points in the two sets */ 1223 static PetscErrorCode PetscQuadraturePointsMerge(PetscQuadrature quadA, PetscQuadrature quadB, PetscQuadrature *quadJoint, PetscInt *aToJoint[], PetscInt *bToJoint[]) 1224 { 1225 PetscInt dimA, dimB; 1226 PetscInt nA, nB, nJoint, i, j, d; 1227 const PetscReal *pointsA; 1228 const PetscReal *pointsB; 1229 PetscReal *pointsJoint; 1230 PetscInt *aToJ, *bToJ; 1231 PetscQuadrature qJ; 1232 PetscErrorCode ierr; 1233 1234 PetscFunctionBegin; 1235 ierr = PetscQuadratureGetData(quadA, &dimA, NULL, &nA, &pointsA, NULL);CHKERRQ(ierr); 1236 ierr = PetscQuadratureGetData(quadB, &dimB, NULL, &nB, &pointsB, NULL);CHKERRQ(ierr); 1237 if (dimA != dimB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Quadrature points must be in the same dimension"); 1238 nJoint = nA; 1239 ierr = PetscMalloc1(nA, &aToJ);CHKERRQ(ierr); 1240 for (i = 0; i < nA; i++) aToJ[i] = i; 1241 ierr = PetscMalloc1(nB, &bToJ);CHKERRQ(ierr); 1242 for (i = 0; i < nB; i++) { 1243 for (j = 0; j < nA; j++) { 1244 bToJ[i] = -1; 1245 for (d = 0; d < dimA; d++) if (PetscAbsReal(pointsB[i * dimA + d] - pointsA[j * dimA + d]) > PETSC_SMALL) break; 1246 if (d == dimA) { 1247 bToJ[i] = j; 1248 break; 1249 } 1250 } 1251 if (bToJ[i] == -1) { 1252 bToJ[i] = nJoint++; 1253 } 1254 } 1255 *aToJoint = aToJ; 1256 *bToJoint = bToJ; 1257 ierr = PetscMalloc1(nJoint * dimA, &pointsJoint);CHKERRQ(ierr); 1258 ierr = PetscArraycpy(pointsJoint, pointsA, nA * dimA);CHKERRQ(ierr); 1259 for (i = 0; i < nB; i++) { 1260 if (bToJ[i] >= nA) { 1261 for (d = 0; d < dimA; d++) pointsJoint[bToJ[i] * dimA + d] = pointsB[i * dimA + d]; 1262 } 1263 } 1264 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &qJ);CHKERRQ(ierr); 1265 ierr = PetscQuadratureSetData(qJ, dimA, 0, nJoint, pointsJoint, NULL);CHKERRQ(ierr); 1266 *quadJoint = qJ; 1267 PetscFunctionReturn(0); 1268 } 1269 1270 /* Matrices matA and matB are both quadrature -> dof matrices: produce a matrix that is joint quadrature -> union of 1271 * dofs, where the joint quadrature was produced by PetscQuadraturePointsMerge */ 1272 static PetscErrorCode MatricesMerge(Mat matA, Mat matB, PetscInt dim, PetscInt k, PetscInt numMerged, const PetscInt aToMerged[], const PetscInt bToMerged[], Mat *matMerged) 1273 { 1274 PetscInt m, n, mA, nA, mB, nB, Nk, i, j, l; 1275 Mat M; 1276 PetscInt *nnz; 1277 PetscInt maxnnz; 1278 PetscInt *work; 1279 PetscErrorCode ierr; 1280 1281 PetscFunctionBegin; 1282 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 1283 ierr = MatGetSize(matA, &mA, &nA);CHKERRQ(ierr); 1284 if (nA % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matA column space not a multiple of k-form size"); 1285 ierr = MatGetSize(matB, &mB, &nB);CHKERRQ(ierr); 1286 if (nB % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matB column space not a multiple of k-form size"); 1287 m = mA + mB; 1288 n = numMerged * Nk; 1289 ierr = PetscMalloc1(m, &nnz);CHKERRQ(ierr); 1290 maxnnz = 0; 1291 for (i = 0; i < mA; i++) { 1292 ierr = MatGetRow(matA, i, &(nnz[i]), NULL, NULL);CHKERRQ(ierr); 1293 if (nnz[i] % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matA are not in k-form size blocks"); 1294 maxnnz = PetscMax(maxnnz, nnz[i]); 1295 } 1296 for (i = 0; i < mB; i++) { 1297 ierr = MatGetRow(matB, i, &(nnz[i+mA]), NULL, NULL);CHKERRQ(ierr); 1298 if (nnz[i+mA] % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matB are not in k-form size blocks"); 1299 maxnnz = PetscMax(maxnnz, nnz[i+mA]); 1300 } 1301 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &M);CHKERRQ(ierr); 1302 ierr = PetscFree(nnz);CHKERRQ(ierr); 1303 /* reasoning about which points each dof needs depends on having zeros computed at points preserved */ 1304 ierr = MatSetOption(M, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr); 1305 ierr = PetscMalloc1(maxnnz, &work);CHKERRQ(ierr); 1306 for (i = 0; i < mA; i++) { 1307 const PetscInt *cols; 1308 const PetscScalar *vals; 1309 PetscInt nCols; 1310 ierr = MatGetRow(matA, i, &nCols, &cols, &vals);CHKERRQ(ierr); 1311 for (j = 0; j < nCols / Nk; j++) { 1312 PetscInt newCol = aToMerged[cols[j * Nk] / Nk]; 1313 for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l; 1314 } 1315 ierr = MatSetValuesBlocked(M, 1, &i, nCols, work, vals, INSERT_VALUES);CHKERRQ(ierr); 1316 ierr = MatRestoreRow(matA, i, &nCols, &cols, &vals);CHKERRQ(ierr); 1317 } 1318 for (i = 0; i < mB; i++) { 1319 const PetscInt *cols; 1320 const PetscScalar *vals; 1321 1322 PetscInt row = i + mA; 1323 PetscInt nCols; 1324 ierr = MatGetRow(matB, i, &nCols, &cols, &vals);CHKERRQ(ierr); 1325 for (j = 0; j < nCols / Nk; j++) { 1326 PetscInt newCol = bToMerged[cols[j * Nk] / Nk]; 1327 for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l; 1328 } 1329 ierr = MatSetValuesBlocked(M, 1, &row, nCols, work, vals, INSERT_VALUES);CHKERRQ(ierr); 1330 ierr = MatRestoreRow(matB, i, &nCols, &cols, &vals);CHKERRQ(ierr); 1331 } 1332 ierr = PetscFree(work);CHKERRQ(ierr); 1333 ierr = MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1334 ierr = MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1335 *matMerged = M; 1336 PetscFunctionReturn(0); 1337 } 1338 1339 /* Take a dual space and product a segment space that has all the same specifications (trimmed, continuous, order, 1340 * node set), except for the form degree. For computing boundary dofs and for making tensor product spaces */ 1341 static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Lagrange(PetscDualSpace sp, DM K, PetscInt f, PetscInt k, PetscInt Ncopies, PetscBool interiorOnly, PetscDualSpace *bdsp) 1342 { 1343 PetscInt Nknew, Ncnew; 1344 PetscInt dim, pointDim = -1; 1345 PetscInt depth; 1346 DM dm; 1347 PetscDualSpace_Lag *newlag; 1348 PetscErrorCode ierr; 1349 1350 PetscFunctionBegin; 1351 ierr = PetscDualSpaceGetDM(sp,&dm);CHKERRQ(ierr); 1352 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 1353 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 1354 ierr = PetscDualSpaceDuplicate(sp,bdsp);CHKERRQ(ierr); 1355 ierr = PetscDualSpaceSetFormDegree(*bdsp,k);CHKERRQ(ierr); 1356 if (!K) { 1357 PetscBool isSimplex; 1358 1359 if (depth == dim) { 1360 PetscInt coneSize; 1361 1362 pointDim = dim - 1; 1363 ierr = DMPlexGetConeSize(dm,f,&coneSize);CHKERRQ(ierr); 1364 isSimplex = (PetscBool) (coneSize == dim); 1365 ierr = PetscDualSpaceCreateReferenceCell(*bdsp, dim-1, isSimplex, &K);CHKERRQ(ierr); 1366 } else if (depth == 1) { 1367 pointDim = 0; 1368 ierr = PetscDualSpaceCreateReferenceCell(*bdsp, 0, PETSC_TRUE, &K);CHKERRQ(ierr); 1369 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element"); 1370 } else { 1371 ierr = PetscObjectReference((PetscObject)K);CHKERRQ(ierr); 1372 ierr = DMGetDimension(K, &pointDim);CHKERRQ(ierr); 1373 } 1374 ierr = PetscDualSpaceSetDM(*bdsp, K);CHKERRQ(ierr); 1375 ierr = DMDestroy(&K);CHKERRQ(ierr); 1376 ierr = PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew);CHKERRQ(ierr); 1377 Ncnew = Nknew * Ncopies; 1378 ierr = PetscDualSpaceSetNumComponents(*bdsp, Ncnew);CHKERRQ(ierr); 1379 newlag = (PetscDualSpace_Lag *) (*bdsp)->data; 1380 newlag->interiorOnly = interiorOnly; 1381 ierr = PetscDualSpaceSetUp(*bdsp);CHKERRQ(ierr); 1382 PetscFunctionReturn(0); 1383 } 1384 1385 /* Construct simplex nodes from a nodefamily, add Nk dof vectors of length Nk at each node. 1386 * Return the (quadrature, matrix) form of the dofs and the nodeIndices form as well. 1387 * 1388 * Sometimes we want a set of nodes to be contained in the interior of the element, 1389 * even when the node scheme puts nodes on the boundaries. numNodeSkip tells 1390 * the routine how many "layers" of nodes need to be skipped. 1391 * */ 1392 static PetscErrorCode PetscDualSpaceLagrangeCreateSimplexNodeMat(Petsc1DNodeFamily nodeFamily, PetscInt dim, PetscInt sum, PetscInt Nk, PetscInt numNodeSkip, PetscQuadrature *iNodes, Mat *iMat, PetscLagNodeIndices *nodeIndices) 1393 { 1394 PetscReal *extraNodeCoords, *nodeCoords; 1395 PetscInt nNodes, nExtraNodes; 1396 PetscInt i, j, k, extraSum = sum + numNodeSkip * (1 + dim); 1397 PetscQuadrature intNodes; 1398 Mat intMat; 1399 PetscLagNodeIndices ni; 1400 PetscErrorCode ierr; 1401 1402 PetscFunctionBegin; 1403 ierr = PetscDTBinomialInt(dim + sum, dim, &nNodes);CHKERRQ(ierr); 1404 ierr = PetscDTBinomialInt(dim + extraSum, dim, &nExtraNodes);CHKERRQ(ierr); 1405 1406 ierr = PetscMalloc1(dim * nExtraNodes, &extraNodeCoords);CHKERRQ(ierr); 1407 ierr = PetscNew(&ni);CHKERRQ(ierr); 1408 ni->nodeIdxDim = dim + 1; 1409 ni->nodeVecDim = Nk; 1410 ni->nNodes = nNodes * Nk; 1411 ni->refct = 1; 1412 ierr = PetscMalloc1(nNodes * Nk * (dim + 1), &(ni->nodeIdx));CHKERRQ(ierr); 1413 ierr = PetscMalloc1(nNodes * Nk * Nk, &(ni->nodeVec));CHKERRQ(ierr); 1414 for (i = 0; i < nNodes; i++) for (j = 0; j < Nk; j++) for (k = 0; k < Nk; k++) ni->nodeVec[(i * Nk + j) * Nk + k] = (j == k) ? 1. : 0.; 1415 ierr = Petsc1DNodeFamilyComputeSimplexNodes(nodeFamily, dim, extraSum, extraNodeCoords);CHKERRQ(ierr); 1416 if (numNodeSkip) { 1417 PetscInt k; 1418 PetscInt *tup; 1419 1420 ierr = PetscMalloc1(dim * nNodes, &nodeCoords);CHKERRQ(ierr); 1421 ierr = PetscMalloc1(dim + 1, &tup);CHKERRQ(ierr); 1422 for (k = 0; k < nNodes; k++) { 1423 PetscInt j, c; 1424 PetscInt index; 1425 1426 ierr = PetscDTIndexToBary(dim + 1, sum, k, tup);CHKERRQ(ierr); 1427 for (j = 0; j < dim + 1; j++) tup[j] += numNodeSkip; 1428 for (c = 0; c < Nk; c++) { 1429 for (j = 0; j < dim + 1; j++) { 1430 ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1; 1431 } 1432 } 1433 ierr = PetscDTBaryToIndex(dim + 1, extraSum, tup, &index);CHKERRQ(ierr); 1434 for (j = 0; j < dim; j++) nodeCoords[k * dim + j] = extraNodeCoords[index * dim + j]; 1435 } 1436 ierr = PetscFree(tup);CHKERRQ(ierr); 1437 ierr = PetscFree(extraNodeCoords);CHKERRQ(ierr); 1438 } else { 1439 PetscInt k; 1440 PetscInt *tup; 1441 1442 nodeCoords = extraNodeCoords; 1443 ierr = PetscMalloc1(dim + 1, &tup);CHKERRQ(ierr); 1444 for (k = 0; k < nNodes; k++) { 1445 PetscInt j, c; 1446 1447 ierr = PetscDTIndexToBary(dim + 1, sum, k, tup);CHKERRQ(ierr); 1448 for (c = 0; c < Nk; c++) { 1449 for (j = 0; j < dim + 1; j++) { 1450 /* barycentric indices can have zeros, but we don't want to push forward zeros because it makes it harder to 1451 * determine which nodes correspond to which under symmetries, so we increase by 1. This is fine 1452 * because the nodeIdx coordinates don't have any meaning other than helping to identify symmetries */ 1453 ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1; 1454 } 1455 } 1456 } 1457 ierr = PetscFree(tup);CHKERRQ(ierr); 1458 } 1459 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &intNodes);CHKERRQ(ierr); 1460 ierr = PetscQuadratureSetData(intNodes, dim, 0, nNodes, nodeCoords, NULL);CHKERRQ(ierr); 1461 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes * Nk, nNodes * Nk, Nk, NULL, &intMat);CHKERRQ(ierr); 1462 ierr = MatSetOption(intMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);CHKERRQ(ierr); 1463 for (j = 0; j < nNodes * Nk; j++) { 1464 PetscInt rem = j % Nk; 1465 PetscInt a, aprev = j - rem; 1466 PetscInt anext = aprev + Nk; 1467 1468 for (a = aprev; a < anext; a++) { 1469 ierr = MatSetValue(intMat, j, a, (a == j) ? 1. : 0., INSERT_VALUES);CHKERRQ(ierr); 1470 } 1471 } 1472 ierr = MatAssemblyBegin(intMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1473 ierr = MatAssemblyEnd(intMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1474 *iNodes = intNodes; 1475 *iMat = intMat; 1476 *nodeIndices = ni; 1477 PetscFunctionReturn(0); 1478 } 1479 1480 /* once the nodeIndices have been created for the interior of the reference cell, and for all of the boundary cells, 1481 * push forward the boundary dofs and concatenate them into the full node indices for the dual space */ 1482 static PetscErrorCode PetscDualSpaceLagrangeCreateAllNodeIdx(PetscDualSpace sp) 1483 { 1484 DM dm; 1485 PetscInt dim, nDofs; 1486 PetscSection section; 1487 PetscInt pStart, pEnd, p; 1488 PetscInt formDegree, Nk; 1489 PetscInt nodeIdxDim, spintdim; 1490 PetscDualSpace_Lag *lag; 1491 PetscLagNodeIndices ni, verti; 1492 PetscErrorCode ierr; 1493 1494 PetscFunctionBegin; 1495 lag = (PetscDualSpace_Lag *) sp->data; 1496 verti = lag->vertIndices; 1497 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 1498 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1499 ierr = PetscDualSpaceGetFormDegree(sp, &formDegree);CHKERRQ(ierr); 1500 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);CHKERRQ(ierr); 1501 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 1502 ierr = PetscSectionGetStorageSize(section, &nDofs);CHKERRQ(ierr); 1503 ierr = PetscNew(&ni);CHKERRQ(ierr); 1504 ni->nodeIdxDim = nodeIdxDim = verti->nodeIdxDim; 1505 ni->nodeVecDim = Nk; 1506 ni->nNodes = nDofs; 1507 ni->refct = 1; 1508 ierr = PetscMalloc1(nodeIdxDim * nDofs, &(ni->nodeIdx));CHKERRQ(ierr); 1509 ierr = PetscMalloc1(Nk * nDofs, &(ni->nodeVec));CHKERRQ(ierr); 1510 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1511 ierr = PetscSectionGetDof(section, 0, &spintdim);CHKERRQ(ierr); 1512 if (spintdim) { 1513 ierr = PetscArraycpy(ni->nodeIdx, lag->intNodeIndices->nodeIdx, spintdim * nodeIdxDim);CHKERRQ(ierr); 1514 ierr = PetscArraycpy(ni->nodeVec, lag->intNodeIndices->nodeVec, spintdim * Nk);CHKERRQ(ierr); 1515 } 1516 for (p = pStart + 1; p < pEnd; p++) { 1517 PetscDualSpace psp = sp->pointSpaces[p]; 1518 PetscDualSpace_Lag *plag; 1519 PetscInt dof, off; 1520 1521 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 1522 if (!dof) continue; 1523 plag = (PetscDualSpace_Lag *) psp->data; 1524 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 1525 ierr = PetscLagNodeIndicesPushForward(dm, verti, p, plag->vertIndices, plag->intNodeIndices, 0, formDegree, &(ni->nodeIdx[off * nodeIdxDim]), &(ni->nodeVec[off * Nk]));CHKERRQ(ierr); 1526 } 1527 lag->allNodeIndices = ni; 1528 PetscFunctionReturn(0); 1529 } 1530 1531 /* once the (quadrature, Matrix) forms of the dofs have been created for the interior of the 1532 * reference cell and for the boundary cells, jk 1533 * push forward the boundary data and concatenate them into the full (quadrature, matrix) data 1534 * for the dual space */ 1535 static PetscErrorCode PetscDualSpaceCreateAllDataFromInteriorData(PetscDualSpace sp) 1536 { 1537 DM dm; 1538 PetscSection section; 1539 PetscInt pStart, pEnd, p, k, Nk, dim, Nc; 1540 PetscInt nNodes; 1541 PetscInt countNodes; 1542 Mat allMat; 1543 PetscQuadrature allNodes; 1544 PetscInt nDofs; 1545 PetscInt maxNzforms, j; 1546 PetscScalar *work; 1547 PetscReal *L, *J, *Jinv, *v0, *pv0; 1548 PetscInt *iwork; 1549 PetscReal *nodes; 1550 PetscErrorCode ierr; 1551 1552 PetscFunctionBegin; 1553 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 1554 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1555 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 1556 ierr = PetscSectionGetStorageSize(section, &nDofs);CHKERRQ(ierr); 1557 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1558 ierr = PetscDualSpaceGetFormDegree(sp, &k);CHKERRQ(ierr); 1559 ierr = PetscDualSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 1560 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 1561 for (p = pStart, nNodes = 0, maxNzforms = 0; p < pEnd; p++) { 1562 PetscDualSpace psp; 1563 DM pdm; 1564 PetscInt pdim, pNk; 1565 PetscQuadrature intNodes; 1566 Mat intMat; 1567 1568 ierr = PetscDualSpaceGetPointSubspace(sp, p, &psp);CHKERRQ(ierr); 1569 if (!psp) continue; 1570 ierr = PetscDualSpaceGetDM(psp, &pdm);CHKERRQ(ierr); 1571 ierr = DMGetDimension(pdm, &pdim);CHKERRQ(ierr); 1572 if (pdim < PetscAbsInt(k)) continue; 1573 ierr = PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);CHKERRQ(ierr); 1574 ierr = PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);CHKERRQ(ierr); 1575 if (intNodes) { 1576 PetscInt nNodesp; 1577 1578 ierr = PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, NULL, NULL);CHKERRQ(ierr); 1579 nNodes += nNodesp; 1580 } 1581 if (intMat) { 1582 PetscInt maxNzsp; 1583 PetscInt maxNzformsp; 1584 1585 ierr = MatSeqAIJGetMaxRowNonzeros(intMat, &maxNzsp);CHKERRQ(ierr); 1586 if (maxNzsp % pNk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms"); 1587 maxNzformsp = maxNzsp / pNk; 1588 maxNzforms = PetscMax(maxNzforms, maxNzformsp); 1589 } 1590 } 1591 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nDofs, nNodes * Nc, maxNzforms * Nk, NULL, &allMat);CHKERRQ(ierr); 1592 ierr = MatSetOption(allMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);CHKERRQ(ierr); 1593 ierr = PetscMalloc7(dim, &v0, dim, &pv0, dim * dim, &J, dim * dim, &Jinv, Nk * Nk, &L, maxNzforms * Nk, &work, maxNzforms * Nk, &iwork);CHKERRQ(ierr); 1594 for (j = 0; j < dim; j++) pv0[j] = -1.; 1595 ierr = PetscMalloc1(dim * nNodes, &nodes);CHKERRQ(ierr); 1596 for (p = pStart, countNodes = 0; p < pEnd; p++) { 1597 PetscDualSpace psp; 1598 PetscQuadrature intNodes; 1599 DM pdm; 1600 PetscInt pdim, pNk; 1601 PetscInt countNodesIn = countNodes; 1602 PetscReal detJ; 1603 Mat intMat; 1604 1605 ierr = PetscDualSpaceGetPointSubspace(sp, p, &psp);CHKERRQ(ierr); 1606 if (!psp) continue; 1607 ierr = PetscDualSpaceGetDM(psp, &pdm);CHKERRQ(ierr); 1608 ierr = DMGetDimension(pdm, &pdim);CHKERRQ(ierr); 1609 if (pdim < PetscAbsInt(k)) continue; 1610 ierr = PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);CHKERRQ(ierr); 1611 if (intNodes == NULL && intMat == NULL) continue; 1612 ierr = PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);CHKERRQ(ierr); 1613 if (p) { 1614 ierr = DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, Jinv, &detJ);CHKERRQ(ierr); 1615 } else { /* identity */ 1616 PetscInt i,j; 1617 1618 for (i = 0; i < dim; i++) for (j = 0; j < dim; j++) J[i * dim + j] = Jinv[i * dim + j] = 0.; 1619 for (i = 0; i < dim; i++) J[i * dim + i] = Jinv[i * dim + i] = 1.; 1620 for (i = 0; i < dim; i++) v0[i] = -1.; 1621 } 1622 if (pdim != dim) { /* compactify Jacobian */ 1623 PetscInt i, j; 1624 1625 for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j]; 1626 } 1627 ierr = PetscDTAltVPullbackMatrix(pdim, dim, J, k, L);CHKERRQ(ierr); 1628 if (intNodes) { /* push forward quadrature locations by the affine transformation */ 1629 PetscInt nNodesp; 1630 const PetscReal *nodesp; 1631 PetscInt j; 1632 1633 ierr = PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, &nodesp, NULL);CHKERRQ(ierr); 1634 for (j = 0; j < nNodesp; j++, countNodes++) { 1635 PetscInt d, e; 1636 1637 for (d = 0; d < dim; d++) { 1638 nodes[countNodes * dim + d] = v0[d]; 1639 for (e = 0; e < pdim; e++) { 1640 nodes[countNodes * dim + d] += J[d * pdim + e] * (nodesp[j * pdim + e] - pv0[e]); 1641 } 1642 } 1643 } 1644 } 1645 if (intMat) { 1646 PetscInt nrows; 1647 PetscInt off; 1648 1649 ierr = PetscSectionGetDof(section, p, &nrows);CHKERRQ(ierr); 1650 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 1651 for (j = 0; j < nrows; j++) { 1652 PetscInt ncols; 1653 const PetscInt *cols; 1654 const PetscScalar *vals; 1655 PetscInt l, d, e; 1656 PetscInt row = j + off; 1657 1658 ierr = MatGetRow(intMat, j, &ncols, &cols, &vals);CHKERRQ(ierr); 1659 if (ncols % pNk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms"); 1660 for (l = 0; l < ncols / pNk; l++) { 1661 PetscInt blockcol; 1662 1663 for (d = 0; d < pNk; d++) { 1664 if ((cols[l * pNk + d] % pNk) != d) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms"); 1665 } 1666 blockcol = cols[l * pNk] / pNk; 1667 for (d = 0; d < Nk; d++) { 1668 iwork[l * Nk + d] = (blockcol + countNodesIn) * Nk + d; 1669 } 1670 for (d = 0; d < Nk; d++) work[l * Nk + d] = 0.; 1671 for (d = 0; d < Nk; d++) { 1672 for (e = 0; e < pNk; e++) { 1673 /* "push forward" dof by pulling back a k-form to be evaluated on the point: multiply on the right by L */ 1674 work[l * Nk + d] += vals[l * pNk + e] * L[e * Nk + d]; 1675 } 1676 } 1677 } 1678 ierr = MatSetValues(allMat, 1, &row, (ncols / pNk) * Nk, iwork, work, INSERT_VALUES);CHKERRQ(ierr); 1679 ierr = MatRestoreRow(intMat, j, &ncols, &cols, &vals);CHKERRQ(ierr); 1680 } 1681 } 1682 } 1683 ierr = MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1684 ierr = MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1685 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &allNodes);CHKERRQ(ierr); 1686 ierr = PetscQuadratureSetData(allNodes, dim, 0, nNodes, nodes, NULL);CHKERRQ(ierr); 1687 ierr = PetscFree7(v0, pv0, J, Jinv, L, work, iwork);CHKERRQ(ierr); 1688 ierr = MatDestroy(&(sp->allMat));CHKERRQ(ierr); 1689 sp->allMat = allMat; 1690 ierr = PetscQuadratureDestroy(&(sp->allNodes));CHKERRQ(ierr); 1691 sp->allNodes = allNodes; 1692 PetscFunctionReturn(0); 1693 } 1694 1695 /* rather than trying to get all data from the functionals, we create 1696 * the functionals from rows of the quadrature -> dof matrix. 1697 * 1698 * Ideally most of the uses of PetscDualSpace in PetscFE will switch 1699 * to using intMat and allMat, so that the individual functionals 1700 * don't need to be constructed at all */ 1701 static PetscErrorCode PetscDualSpaceComputeFunctionalsFromAllData(PetscDualSpace sp) 1702 { 1703 PetscQuadrature allNodes; 1704 Mat allMat; 1705 PetscInt nDofs; 1706 PetscInt dim, k, Nk, Nc, f; 1707 DM dm; 1708 PetscInt nNodes, spdim; 1709 const PetscReal *nodes = NULL; 1710 PetscSection section; 1711 PetscBool useMoments; 1712 PetscErrorCode ierr; 1713 1714 PetscFunctionBegin; 1715 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 1716 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1717 ierr = PetscDualSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 1718 ierr = PetscDualSpaceGetFormDegree(sp, &k);CHKERRQ(ierr); 1719 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 1720 ierr = PetscDualSpaceGetAllData(sp, &allNodes, &allMat);CHKERRQ(ierr); 1721 nNodes = 0; 1722 if (allNodes) { 1723 ierr = PetscQuadratureGetData(allNodes, NULL, NULL, &nNodes, &nodes, NULL);CHKERRQ(ierr); 1724 } 1725 ierr = MatGetSize(allMat, &nDofs, NULL);CHKERRQ(ierr); 1726 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 1727 ierr = PetscSectionGetStorageSize(section, &spdim);CHKERRQ(ierr); 1728 if (spdim != nDofs) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "incompatible all matrix size"); 1729 ierr = PetscMalloc1(nDofs, &(sp->functional));CHKERRQ(ierr); 1730 ierr = PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments);CHKERRQ(ierr); 1731 if (useMoments) { 1732 Mat allMat; 1733 PetscInt momentOrder, i; 1734 PetscBool tensor; 1735 const PetscReal *weights; 1736 PetscScalar *array; 1737 1738 if (nDofs != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "We do not yet support moments beyond P0, nDofs == %D", nDofs); 1739 ierr = PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder);CHKERRQ(ierr); 1740 ierr = PetscDualSpaceLagrangeGetTensor(sp, &tensor);CHKERRQ(ierr); 1741 if (!tensor) {ierr = PetscDTStroudConicalQuadrature(dim, Nc, PetscMax(momentOrder + 1,1), -1.0, 1.0, &(sp->functional[0]));CHKERRQ(ierr);} 1742 else {ierr = PetscDTGaussTensorQuadrature(dim, Nc, PetscMax(momentOrder + 1,1), -1.0, 1.0, &(sp->functional[0]));CHKERRQ(ierr);} 1743 /* Need to replace allNodes and allMat */ 1744 ierr = PetscObjectReference((PetscObject) sp->functional[0]);CHKERRQ(ierr); 1745 ierr = PetscQuadratureDestroy(&(sp->allNodes));CHKERRQ(ierr); 1746 sp->allNodes = sp->functional[0]; 1747 ierr = PetscQuadratureGetData(sp->allNodes, NULL, NULL, &nNodes, NULL, &weights);CHKERRQ(ierr); 1748 ierr = MatCreateSeqDense(PETSC_COMM_SELF, nDofs, nNodes * Nc, NULL, &allMat);CHKERRQ(ierr); 1749 ierr = MatDenseGetArrayWrite(allMat, &array);CHKERRQ(ierr); 1750 for (i = 0; i < nNodes * Nc; ++i) array[i] = weights[i]; 1751 ierr = MatDenseRestoreArrayWrite(allMat, &array);CHKERRQ(ierr); 1752 ierr = MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1753 ierr = MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1754 ierr = MatDestroy(&(sp->allMat));CHKERRQ(ierr); 1755 sp->allMat = allMat; 1756 PetscFunctionReturn(0); 1757 } 1758 for (f = 0; f < nDofs; f++) { 1759 PetscInt ncols, c; 1760 const PetscInt *cols; 1761 const PetscScalar *vals; 1762 PetscReal *nodesf; 1763 PetscReal *weightsf; 1764 PetscInt nNodesf; 1765 PetscInt countNodes; 1766 1767 ierr = MatGetRow(allMat, f, &ncols, &cols, &vals);CHKERRQ(ierr); 1768 if (ncols % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "all matrix is not laid out as blocks of k-forms"); 1769 for (c = 1, nNodesf = 1; c < ncols; c++) { 1770 if ((cols[c] / Nc) != (cols[c-1] / Nc)) nNodesf++; 1771 } 1772 ierr = PetscMalloc1(dim * nNodesf, &nodesf);CHKERRQ(ierr); 1773 ierr = PetscMalloc1(Nc * nNodesf, &weightsf);CHKERRQ(ierr); 1774 for (c = 0, countNodes = 0; c < ncols; c++) { 1775 if (!c || ((cols[c] / Nc) != (cols[c-1] / Nc))) { 1776 PetscInt d; 1777 1778 for (d = 0; d < Nc; d++) { 1779 weightsf[countNodes * Nc + d] = 0.; 1780 } 1781 for (d = 0; d < dim; d++) { 1782 nodesf[countNodes * dim + d] = nodes[(cols[c] / Nc) * dim + d]; 1783 } 1784 countNodes++; 1785 } 1786 weightsf[(countNodes - 1) * Nc + (cols[c] % Nc)] = PetscRealPart(vals[c]); 1787 } 1788 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &(sp->functional[f]));CHKERRQ(ierr); 1789 ierr = PetscQuadratureSetData(sp->functional[f], dim, Nc, nNodesf, nodesf, weightsf);CHKERRQ(ierr); 1790 ierr = MatRestoreRow(allMat, f, &ncols, &cols, &vals);CHKERRQ(ierr); 1791 } 1792 PetscFunctionReturn(0); 1793 } 1794 1795 /* take a matrix meant for k-forms and expand it to one for Ncopies */ 1796 static PetscErrorCode PetscDualSpaceLagrangeMatrixCreateCopies(Mat A, PetscInt Nk, PetscInt Ncopies, Mat *Abs) 1797 { 1798 PetscInt m, n, i, j, k; 1799 PetscInt maxnnz, *nnz, *iwork; 1800 Mat Ac; 1801 PetscErrorCode ierr; 1802 1803 PetscFunctionBegin; 1804 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 1805 if (n % Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of columns in A %D is not a multiple of Nk %D", n, Nk); 1806 ierr = PetscMalloc1(m * Ncopies, &nnz);CHKERRQ(ierr); 1807 for (i = 0, maxnnz = 0; i < m; i++) { 1808 PetscInt innz; 1809 ierr = MatGetRow(A, i, &innz, NULL, NULL);CHKERRQ(ierr); 1810 if (innz % Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "A row %D nnzs is not a multiple of Nk %D", innz, Nk); 1811 for (j = 0; j < Ncopies; j++) nnz[i * Ncopies + j] = innz; 1812 maxnnz = PetscMax(maxnnz, innz); 1813 } 1814 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, m * Ncopies, n * Ncopies, 0, nnz, &Ac);CHKERRQ(ierr); 1815 ierr = MatSetOption(Ac, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr); 1816 ierr = PetscFree(nnz);CHKERRQ(ierr); 1817 ierr = PetscMalloc1(maxnnz, &iwork);CHKERRQ(ierr); 1818 for (i = 0; i < m; i++) { 1819 PetscInt innz; 1820 const PetscInt *cols; 1821 const PetscScalar *vals; 1822 1823 ierr = MatGetRow(A, i, &innz, &cols, &vals);CHKERRQ(ierr); 1824 for (j = 0; j < innz; j++) iwork[j] = (cols[j] / Nk) * (Nk * Ncopies) + (cols[j] % Nk); 1825 for (j = 0; j < Ncopies; j++) { 1826 PetscInt row = i * Ncopies + j; 1827 1828 ierr = MatSetValues(Ac, 1, &row, innz, iwork, vals, INSERT_VALUES);CHKERRQ(ierr); 1829 for (k = 0; k < innz; k++) iwork[k] += Nk; 1830 } 1831 ierr = MatRestoreRow(A, i, &innz, &cols, &vals);CHKERRQ(ierr); 1832 } 1833 ierr = PetscFree(iwork);CHKERRQ(ierr); 1834 ierr = MatAssemblyBegin(Ac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1835 ierr = MatAssemblyEnd(Ac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1836 *Abs = Ac; 1837 PetscFunctionReturn(0); 1838 } 1839 1840 /* check if a cell is a tensor product of the segment with a facet, 1841 * specifically checking if f and f2 can be the "endpoints" (like the triangles 1842 * at either end of a wedge) */ 1843 static PetscErrorCode DMPlexPointIsTensor_Internal_Given(DM dm, PetscInt p, PetscInt f, PetscInt f2, PetscBool *isTensor) 1844 { 1845 PetscInt coneSize, c; 1846 const PetscInt *cone; 1847 const PetscInt *fCone; 1848 const PetscInt *f2Cone; 1849 PetscInt fs[2]; 1850 PetscInt meetSize, nmeet; 1851 const PetscInt *meet; 1852 PetscErrorCode ierr; 1853 1854 PetscFunctionBegin; 1855 fs[0] = f; 1856 fs[1] = f2; 1857 ierr = DMPlexGetMeet(dm, 2, fs, &meetSize, &meet);CHKERRQ(ierr); 1858 nmeet = meetSize; 1859 ierr = DMPlexRestoreMeet(dm, 2, fs, &meetSize, &meet);CHKERRQ(ierr); 1860 /* two points that have a non-empty meet cannot be at opposite ends of a cell */ 1861 if (nmeet) { 1862 *isTensor = PETSC_FALSE; 1863 PetscFunctionReturn(0); 1864 } 1865 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1866 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 1867 ierr = DMPlexGetCone(dm, f, &fCone);CHKERRQ(ierr); 1868 ierr = DMPlexGetCone(dm, f2, &f2Cone);CHKERRQ(ierr); 1869 for (c = 0; c < coneSize; c++) { 1870 PetscInt e, ef; 1871 PetscInt d = -1, d2 = -1; 1872 PetscInt dcount, d2count; 1873 PetscInt t = cone[c]; 1874 PetscInt tConeSize; 1875 PetscBool tIsTensor; 1876 const PetscInt *tCone; 1877 1878 if (t == f || t == f2) continue; 1879 /* for every other facet in the cone, check that is has 1880 * one ridge in common with each end */ 1881 ierr = DMPlexGetConeSize(dm, t, &tConeSize);CHKERRQ(ierr); 1882 ierr = DMPlexGetCone(dm, t, &tCone);CHKERRQ(ierr); 1883 1884 dcount = 0; 1885 d2count = 0; 1886 for (e = 0; e < tConeSize; e++) { 1887 PetscInt q = tCone[e]; 1888 for (ef = 0; ef < coneSize - 2; ef++) { 1889 if (fCone[ef] == q) { 1890 if (dcount) { 1891 *isTensor = PETSC_FALSE; 1892 PetscFunctionReturn(0); 1893 } 1894 d = q; 1895 dcount++; 1896 } else if (f2Cone[ef] == q) { 1897 if (d2count) { 1898 *isTensor = PETSC_FALSE; 1899 PetscFunctionReturn(0); 1900 } 1901 d2 = q; 1902 d2count++; 1903 } 1904 } 1905 } 1906 /* if the whole cell is a tensor with the segment, then this 1907 * facet should be a tensor with the segment */ 1908 ierr = DMPlexPointIsTensor_Internal_Given(dm, t, d, d2, &tIsTensor);CHKERRQ(ierr); 1909 if (!tIsTensor) { 1910 *isTensor = PETSC_FALSE; 1911 PetscFunctionReturn(0); 1912 } 1913 } 1914 *isTensor = PETSC_TRUE; 1915 PetscFunctionReturn(0); 1916 } 1917 1918 /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair 1919 * that could be the opposite ends */ 1920 static PetscErrorCode DMPlexPointIsTensor_Internal(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB) 1921 { 1922 PetscInt coneSize, c, c2; 1923 const PetscInt *cone; 1924 PetscErrorCode ierr; 1925 1926 PetscFunctionBegin; 1927 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1928 if (!coneSize) { 1929 if (isTensor) *isTensor = PETSC_FALSE; 1930 if (endA) *endA = -1; 1931 if (endB) *endB = -1; 1932 } 1933 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 1934 for (c = 0; c < coneSize; c++) { 1935 PetscInt f = cone[c]; 1936 PetscInt fConeSize; 1937 1938 ierr = DMPlexGetConeSize(dm, f, &fConeSize);CHKERRQ(ierr); 1939 if (fConeSize != coneSize - 2) continue; 1940 1941 for (c2 = c + 1; c2 < coneSize; c2++) { 1942 PetscInt f2 = cone[c2]; 1943 PetscBool isTensorff2; 1944 PetscInt f2ConeSize; 1945 1946 ierr = DMPlexGetConeSize(dm, f2, &f2ConeSize);CHKERRQ(ierr); 1947 if (f2ConeSize != coneSize - 2) continue; 1948 1949 ierr = DMPlexPointIsTensor_Internal_Given(dm, p, f, f2, &isTensorff2);CHKERRQ(ierr); 1950 if (isTensorff2) { 1951 if (isTensor) *isTensor = PETSC_TRUE; 1952 if (endA) *endA = f; 1953 if (endB) *endB = f2; 1954 PetscFunctionReturn(0); 1955 } 1956 } 1957 } 1958 if (isTensor) *isTensor = PETSC_FALSE; 1959 if (endA) *endA = -1; 1960 if (endB) *endB = -1; 1961 PetscFunctionReturn(0); 1962 } 1963 1964 /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair 1965 * that could be the opposite ends */ 1966 static PetscErrorCode DMPlexPointIsTensor(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB) 1967 { 1968 DMPlexInterpolatedFlag interpolated; 1969 PetscErrorCode ierr; 1970 1971 PetscFunctionBegin; 1972 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1973 if (interpolated != DMPLEX_INTERPOLATED_FULL) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Only for interpolated DMPlex's"); 1974 ierr = DMPlexPointIsTensor_Internal(dm, p, isTensor, endA, endB);CHKERRQ(ierr); 1975 PetscFunctionReturn(0); 1976 } 1977 1978 /* Let k = formDegree and k' = -sign(k) * dim + k. Transform a symmetric frame for k-forms on the biunit simplex into 1979 * a symmetric frame for k'-forms on the biunit simplex. 1980 * 1981 * A frame is "symmetric" if the pullback of every symmetry of the biunit simplex is a permutation of the frame. 1982 * 1983 * forms in the symmetric frame are used as dofs in the untrimmed simplex spaces. This way, symmetries of the 1984 * reference cell result in permutations of dofs grouped by node. 1985 * 1986 * Use T to transform dof matrices for k'-forms into dof matrices for k-forms as a block diagonal transformation on 1987 * the right. 1988 */ 1989 static PetscErrorCode BiunitSimplexSymmetricFormTransformation(PetscInt dim, PetscInt formDegree, PetscReal T[]) 1990 { 1991 PetscInt k = formDegree; 1992 PetscInt kd = k < 0 ? dim + k : k - dim; 1993 PetscInt Nk; 1994 PetscReal *biToEq, *eqToBi, *biToEqStar, *eqToBiStar; 1995 PetscInt fact; 1996 PetscErrorCode ierr; 1997 1998 PetscFunctionBegin; 1999 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 2000 ierr = PetscCalloc4(dim * dim, &biToEq, dim * dim, &eqToBi, Nk * Nk, &biToEqStar, Nk * Nk, &eqToBiStar);CHKERRQ(ierr); 2001 /* fill in biToEq: Jacobian of the transformation from the biunit simplex to the equilateral simplex */ 2002 fact = 0; 2003 for (PetscInt i = 0; i < dim; i++) { 2004 biToEq[i * dim + i] = PetscSqrtReal(((PetscReal)i + 2.) / (2.*((PetscReal)i+1.))); 2005 fact += 4*(i+1); 2006 for (PetscInt j = i+1; j < dim; j++) { 2007 biToEq[i * dim + j] = PetscSqrtReal(1./(PetscReal)fact); 2008 } 2009 } 2010 /* fill in eqToBi: Jacobian of the transformation from the equilateral simplex to the biunit simplex */ 2011 fact = 0; 2012 for (PetscInt j = 0; j < dim; j++) { 2013 eqToBi[j * dim + j] = PetscSqrtReal(2.*((PetscReal)j+1.)/((PetscReal)j+2)); 2014 fact += j+1; 2015 for (PetscInt i = 0; i < j; i++) { 2016 eqToBi[i * dim + j] = -PetscSqrtReal(1./(PetscReal)fact); 2017 } 2018 } 2019 ierr = PetscDTAltVPullbackMatrix(dim, dim, biToEq, kd, biToEqStar);CHKERRQ(ierr); 2020 ierr = PetscDTAltVPullbackMatrix(dim, dim, eqToBi, k, eqToBiStar);CHKERRQ(ierr); 2021 /* product of pullbacks simulates the following steps 2022 * 2023 * 1. start with frame W = [w_1, w_2, ..., w_m] of k forms that is symmetric on the biunit simplex: 2024 if J is the Jacobian of a symmetry of the biunit simplex, then J_k* W = [J_k*w_1, ..., J_k*w_m] 2025 is a permutation of W. 2026 Even though a k' form --- a (dim - k) form represented by its Hodge star --- has the same geometric 2027 content as a k form, W is not a symmetric frame of k' forms on the biunit simplex. That's because, 2028 for general Jacobian J, J_k* != J_k'*. 2029 * 2. pullback W to the equilateral triangle using the k pullback, W_eq = eqToBi_k* W. All symmetries of the 2030 equilateral simplex have orthonormal Jacobians. For an orthonormal Jacobian O, J_k* = J_k'*, so W_eq is 2031 also a symmetric frame for k' forms on the equilateral simplex. 2032 3. pullback W_eq back to the biunit simplex using the k' pulback, V = biToEq_k'* W_eq = biToEq_k'* eqToBi_k* W. 2033 V is a symmetric frame for k' forms on the biunit simplex. 2034 */ 2035 for (PetscInt i = 0; i < Nk; i++) { 2036 for (PetscInt j = 0; j < Nk; j++) { 2037 PetscReal val = 0.; 2038 for (PetscInt k = 0; k < Nk; k++) val += biToEqStar[i * Nk + k] * eqToBiStar[k * Nk + j]; 2039 T[i * Nk + j] = val; 2040 } 2041 } 2042 ierr = PetscFree4(biToEq, eqToBi, biToEqStar, eqToBiStar);CHKERRQ(ierr); 2043 PetscFunctionReturn(0); 2044 } 2045 2046 /* permute a quadrature -> dof matrix so that its rows are in revlex order by nodeIdx */ 2047 static PetscErrorCode MatPermuteByNodeIdx(Mat A, PetscLagNodeIndices ni, Mat *Aperm) 2048 { 2049 PetscInt m, n, i, j; 2050 PetscInt nodeIdxDim = ni->nodeIdxDim; 2051 PetscInt nodeVecDim = ni->nodeVecDim; 2052 PetscInt *perm; 2053 IS permIS; 2054 IS id; 2055 PetscInt *nIdxPerm; 2056 PetscReal *nVecPerm; 2057 PetscErrorCode ierr; 2058 2059 PetscFunctionBegin; 2060 ierr = PetscLagNodeIndicesGetPermutation(ni, &perm);CHKERRQ(ierr); 2061 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 2062 ierr = PetscMalloc1(nodeIdxDim * m, &nIdxPerm);CHKERRQ(ierr); 2063 ierr = PetscMalloc1(nodeVecDim * m, &nVecPerm);CHKERRQ(ierr); 2064 for (i = 0; i < m; i++) for (j = 0; j < nodeIdxDim; j++) nIdxPerm[i * nodeIdxDim + j] = ni->nodeIdx[perm[i] * nodeIdxDim + j]; 2065 for (i = 0; i < m; i++) for (j = 0; j < nodeVecDim; j++) nVecPerm[i * nodeVecDim + j] = ni->nodeVec[perm[i] * nodeVecDim + j]; 2066 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, perm, PETSC_USE_POINTER, &permIS);CHKERRQ(ierr); 2067 ierr = ISSetPermutation(permIS);CHKERRQ(ierr); 2068 ierr = ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &id);CHKERRQ(ierr); 2069 ierr = ISSetPermutation(id);CHKERRQ(ierr); 2070 ierr = MatPermute(A, permIS, id, Aperm);CHKERRQ(ierr); 2071 ierr = ISDestroy(&permIS);CHKERRQ(ierr); 2072 ierr = ISDestroy(&id);CHKERRQ(ierr); 2073 for (i = 0; i < m; i++) perm[i] = i; 2074 ierr = PetscFree(ni->nodeIdx);CHKERRQ(ierr); 2075 ierr = PetscFree(ni->nodeVec);CHKERRQ(ierr); 2076 ni->nodeIdx = nIdxPerm; 2077 ni->nodeVec = nVecPerm; 2078 PetscFunctionReturn(0); 2079 } 2080 2081 static PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp) 2082 { 2083 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 2084 DM dm = sp->dm; 2085 DM dmint = NULL; 2086 PetscInt order; 2087 PetscInt Nc = sp->Nc; 2088 MPI_Comm comm; 2089 PetscBool continuous; 2090 PetscSection section; 2091 PetscInt depth, dim, pStart, pEnd, cStart, cEnd, p, *pStratStart, *pStratEnd, d; 2092 PetscInt formDegree, Nk, Ncopies; 2093 PetscInt tensorf = -1, tensorf2 = -1; 2094 PetscBool tensorCell, tensorSpace; 2095 PetscBool uniform, trimmed; 2096 Petsc1DNodeFamily nodeFamily; 2097 PetscInt numNodeSkip; 2098 DMPlexInterpolatedFlag interpolated; 2099 PetscBool isbdm; 2100 PetscErrorCode ierr; 2101 2102 PetscFunctionBegin; 2103 /* step 1: sanitize input */ 2104 ierr = PetscObjectGetComm((PetscObject) sp, &comm);CHKERRQ(ierr); 2105 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2106 ierr = PetscObjectTypeCompare((PetscObject)sp, PETSCDUALSPACEBDM, &isbdm);CHKERRQ(ierr); 2107 if (isbdm) { 2108 sp->k = -(dim-1); /* form degree of H-div */ 2109 ierr = PetscObjectChangeTypeName((PetscObject)sp, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 2110 } 2111 ierr = PetscDualSpaceGetFormDegree(sp, &formDegree);CHKERRQ(ierr); 2112 if (PetscAbsInt(formDegree) > dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Form degree must be bounded by dimension"); 2113 ierr = PetscDTBinomialInt(dim,PetscAbsInt(formDegree),&Nk);CHKERRQ(ierr); 2114 if (sp->Nc <= 0 && lag->numCopies > 0) sp->Nc = Nk * lag->numCopies; 2115 Nc = sp->Nc; 2116 if (Nc % Nk) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "Number of components is not a multiple of form degree size"); 2117 if (lag->numCopies <= 0) lag->numCopies = Nc / Nk; 2118 Ncopies = lag->numCopies; 2119 if (Nc / Nk != Ncopies) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "Number of copies * (dim choose k) != Nc"); 2120 if (!dim) sp->order = 0; 2121 order = sp->order; 2122 uniform = sp->uniform; 2123 if (!uniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Variable order not supported yet"); 2124 if (lag->trimmed && !formDegree) lag->trimmed = PETSC_FALSE; /* trimmed spaces are the same as full spaces for 0-forms */ 2125 if (lag->nodeType == PETSCDTNODES_DEFAULT) { 2126 lag->nodeType = PETSCDTNODES_GAUSSJACOBI; 2127 lag->nodeExponent = 0.; 2128 /* trimmed spaces don't include corner vertices, so don't use end nodes by default */ 2129 lag->endNodes = lag->trimmed ? PETSC_FALSE : PETSC_TRUE; 2130 } 2131 /* If a trimmed space and the user did choose nodes with endpoints, skip them by default */ 2132 if (lag->numNodeSkip < 0) lag->numNodeSkip = (lag->trimmed && lag->endNodes) ? 1 : 0; 2133 numNodeSkip = lag->numNodeSkip; 2134 if (lag->trimmed && !order) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zeroth order trimmed elements"); 2135 if (lag->trimmed && PetscAbsInt(formDegree) == dim) { /* convert trimmed n-forms to untrimmed of one polynomial order less */ 2136 sp->order--; 2137 order--; 2138 lag->trimmed = PETSC_FALSE; 2139 } 2140 trimmed = lag->trimmed; 2141 if (!order || PetscAbsInt(formDegree) == dim) lag->continuous = PETSC_FALSE; 2142 continuous = lag->continuous; 2143 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2144 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2145 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2146 if (pStart != 0 || cStart != 0) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Expect DM with chart starting at zero and cells first"); 2147 if (cEnd != 1) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Use PETSCDUALSPACEREFINED for multi-cell reference meshes"); 2148 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 2149 if (interpolated != DMPLEX_INTERPOLATED_FULL) { 2150 ierr = DMPlexInterpolate(dm, &dmint);CHKERRQ(ierr); 2151 } else { 2152 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 2153 dmint = dm; 2154 } 2155 tensorCell = PETSC_FALSE; 2156 if (dim > 1) { 2157 ierr = DMPlexPointIsTensor(dmint, 0, &tensorCell, &tensorf, &tensorf2);CHKERRQ(ierr); 2158 } 2159 lag->tensorCell = tensorCell; 2160 if (dim < 2 || !lag->tensorCell) lag->tensorSpace = PETSC_FALSE; 2161 tensorSpace = lag->tensorSpace; 2162 if (!lag->nodeFamily) { 2163 ierr = Petsc1DNodeFamilyCreate(lag->nodeType, lag->nodeExponent, lag->endNodes, &lag->nodeFamily);CHKERRQ(ierr); 2164 } 2165 nodeFamily = lag->nodeFamily; 2166 if (interpolated != DMPLEX_INTERPOLATED_FULL && continuous && (PetscAbsInt(formDegree) > 0 || order > 1)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Reference element won't support all boundary nodes"); 2167 2168 /* step 2: construct the boundary spaces */ 2169 ierr = PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);CHKERRQ(ierr); 2170 ierr = PetscCalloc1(pEnd,&(sp->pointSpaces));CHKERRQ(ierr); 2171 for (d = 0; d <= depth; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);CHKERRQ(ierr);} 2172 ierr = PetscDualSpaceSectionCreate_Internal(sp, §ion);CHKERRQ(ierr); 2173 sp->pointSection = section; 2174 if (continuous && !(lag->interiorOnly)) { 2175 PetscInt h; 2176 2177 for (p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */ 2178 PetscReal v0[3]; 2179 DMPolytopeType ptype; 2180 PetscReal J[9], detJ; 2181 PetscInt q; 2182 2183 ierr = DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, NULL, &detJ);CHKERRQ(ierr); 2184 ierr = DMPlexGetCellType(dm, p, &ptype);CHKERRQ(ierr); 2185 2186 /* compare to previous facets: if computed, reference that dualspace */ 2187 for (q = pStratStart[depth - 1]; q < p; q++) { 2188 DMPolytopeType qtype; 2189 2190 ierr = DMPlexGetCellType(dm, q, &qtype);CHKERRQ(ierr); 2191 if (qtype == ptype) break; 2192 } 2193 if (q < p) { /* this facet has the same dual space as that one */ 2194 ierr = PetscObjectReference((PetscObject)sp->pointSpaces[q]);CHKERRQ(ierr); 2195 sp->pointSpaces[p] = sp->pointSpaces[q]; 2196 continue; 2197 } 2198 /* if not, recursively compute this dual space */ 2199 ierr = PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,p,formDegree,Ncopies,PETSC_FALSE,&sp->pointSpaces[p]);CHKERRQ(ierr); 2200 } 2201 for (h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */ 2202 PetscInt hd = depth - h; 2203 PetscInt hdim = dim - h; 2204 2205 if (hdim < PetscAbsInt(formDegree)) break; 2206 for (p = pStratStart[hd]; p < pStratEnd[hd]; p++) { 2207 PetscInt suppSize, s; 2208 const PetscInt *supp; 2209 2210 ierr = DMPlexGetSupportSize(dm, p, &suppSize);CHKERRQ(ierr); 2211 ierr = DMPlexGetSupport(dm, p, &supp);CHKERRQ(ierr); 2212 for (s = 0; s < suppSize; s++) { 2213 DM qdm; 2214 PetscDualSpace qsp, psp; 2215 PetscInt c, coneSize, q; 2216 const PetscInt *cone; 2217 const PetscInt *refCone; 2218 2219 q = supp[0]; 2220 qsp = sp->pointSpaces[q]; 2221 ierr = DMPlexGetConeSize(dm, q, &coneSize);CHKERRQ(ierr); 2222 ierr = DMPlexGetCone(dm, q, &cone);CHKERRQ(ierr); 2223 for (c = 0; c < coneSize; c++) if (cone[c] == p) break; 2224 if (c == coneSize) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "cone/support mismatch"); 2225 ierr = PetscDualSpaceGetDM(qsp, &qdm);CHKERRQ(ierr); 2226 ierr = DMPlexGetCone(qdm, 0, &refCone);CHKERRQ(ierr); 2227 /* get the equivalent dual space from the support dual space */ 2228 ierr = PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp);CHKERRQ(ierr); 2229 if (!s) { 2230 ierr = PetscObjectReference((PetscObject)psp);CHKERRQ(ierr); 2231 sp->pointSpaces[p] = psp; 2232 } 2233 } 2234 } 2235 } 2236 for (p = 1; p < pEnd; p++) { 2237 PetscInt pspdim; 2238 if (!sp->pointSpaces[p]) continue; 2239 ierr = PetscDualSpaceGetInteriorDimension(sp->pointSpaces[p], &pspdim);CHKERRQ(ierr); 2240 ierr = PetscSectionSetDof(section, p, pspdim);CHKERRQ(ierr); 2241 } 2242 } 2243 2244 if (Ncopies > 1) { 2245 Mat intMatScalar, allMatScalar; 2246 PetscDualSpace scalarsp; 2247 PetscDualSpace_Lag *scalarlag; 2248 2249 ierr = PetscDualSpaceDuplicate(sp, &scalarsp);CHKERRQ(ierr); 2250 /* Setting the number of components to Nk is a space with 1 copy of each k-form */ 2251 ierr = PetscDualSpaceSetNumComponents(scalarsp, Nk);CHKERRQ(ierr); 2252 ierr = PetscDualSpaceSetUp(scalarsp);CHKERRQ(ierr); 2253 ierr = PetscDualSpaceGetInteriorData(scalarsp, &(sp->intNodes), &intMatScalar);CHKERRQ(ierr); 2254 ierr = PetscObjectReference((PetscObject)(sp->intNodes));CHKERRQ(ierr); 2255 if (intMatScalar) {ierr = PetscDualSpaceLagrangeMatrixCreateCopies(intMatScalar, Nk, Ncopies, &(sp->intMat));CHKERRQ(ierr);} 2256 ierr = PetscDualSpaceGetAllData(scalarsp, &(sp->allNodes), &allMatScalar);CHKERRQ(ierr); 2257 ierr = PetscObjectReference((PetscObject)(sp->allNodes));CHKERRQ(ierr); 2258 ierr = PetscDualSpaceLagrangeMatrixCreateCopies(allMatScalar, Nk, Ncopies, &(sp->allMat));CHKERRQ(ierr); 2259 sp->spdim = scalarsp->spdim * Ncopies; 2260 sp->spintdim = scalarsp->spintdim * Ncopies; 2261 scalarlag = (PetscDualSpace_Lag *) scalarsp->data; 2262 ierr = PetscLagNodeIndicesReference(scalarlag->vertIndices);CHKERRQ(ierr); 2263 lag->vertIndices = scalarlag->vertIndices; 2264 ierr = PetscLagNodeIndicesReference(scalarlag->intNodeIndices);CHKERRQ(ierr); 2265 lag->intNodeIndices = scalarlag->intNodeIndices; 2266 ierr = PetscLagNodeIndicesReference(scalarlag->allNodeIndices);CHKERRQ(ierr); 2267 lag->allNodeIndices = scalarlag->allNodeIndices; 2268 ierr = PetscDualSpaceDestroy(&scalarsp);CHKERRQ(ierr); 2269 ierr = PetscSectionSetDof(section, 0, sp->spintdim);CHKERRQ(ierr); 2270 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 2271 ierr = PetscDualSpaceComputeFunctionalsFromAllData(sp);CHKERRQ(ierr); 2272 ierr = PetscFree2(pStratStart, pStratEnd);CHKERRQ(ierr); 2273 ierr = DMDestroy(&dmint);CHKERRQ(ierr); 2274 PetscFunctionReturn(0); 2275 } 2276 2277 if (trimmed && !continuous) { 2278 /* the dofs of a trimmed space don't have a nice tensor/lattice structure: 2279 * just construct the continuous dual space and copy all of the data over, 2280 * allocating it all to the cell instead of splitting it up between the boundaries */ 2281 PetscDualSpace spcont; 2282 PetscInt spdim, f; 2283 PetscQuadrature allNodes; 2284 PetscDualSpace_Lag *lagc; 2285 Mat allMat; 2286 2287 ierr = PetscDualSpaceDuplicate(sp, &spcont);CHKERRQ(ierr); 2288 ierr = PetscDualSpaceLagrangeSetContinuity(spcont, PETSC_TRUE);CHKERRQ(ierr); 2289 ierr = PetscDualSpaceSetUp(spcont);CHKERRQ(ierr); 2290 ierr = PetscDualSpaceGetDimension(spcont, &spdim);CHKERRQ(ierr); 2291 sp->spdim = sp->spintdim = spdim; 2292 ierr = PetscSectionSetDof(section, 0, spdim);CHKERRQ(ierr); 2293 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 2294 ierr = PetscMalloc1(spdim, &(sp->functional));CHKERRQ(ierr); 2295 for (f = 0; f < spdim; f++) { 2296 PetscQuadrature fn; 2297 2298 ierr = PetscDualSpaceGetFunctional(spcont, f, &fn);CHKERRQ(ierr); 2299 ierr = PetscObjectReference((PetscObject)fn);CHKERRQ(ierr); 2300 sp->functional[f] = fn; 2301 } 2302 ierr = PetscDualSpaceGetAllData(spcont, &allNodes, &allMat);CHKERRQ(ierr); 2303 ierr = PetscObjectReference((PetscObject) allNodes);CHKERRQ(ierr); 2304 ierr = PetscObjectReference((PetscObject) allNodes);CHKERRQ(ierr); 2305 sp->allNodes = sp->intNodes = allNodes; 2306 ierr = PetscObjectReference((PetscObject) allMat);CHKERRQ(ierr); 2307 ierr = PetscObjectReference((PetscObject) allMat);CHKERRQ(ierr); 2308 sp->allMat = sp->intMat = allMat; 2309 lagc = (PetscDualSpace_Lag *) spcont->data; 2310 ierr = PetscLagNodeIndicesReference(lagc->vertIndices);CHKERRQ(ierr); 2311 lag->vertIndices = lagc->vertIndices; 2312 ierr = PetscLagNodeIndicesReference(lagc->allNodeIndices);CHKERRQ(ierr); 2313 ierr = PetscLagNodeIndicesReference(lagc->allNodeIndices);CHKERRQ(ierr); 2314 lag->intNodeIndices = lagc->allNodeIndices; 2315 lag->allNodeIndices = lagc->allNodeIndices; 2316 ierr = PetscDualSpaceDestroy(&spcont);CHKERRQ(ierr); 2317 ierr = PetscFree2(pStratStart, pStratEnd);CHKERRQ(ierr); 2318 ierr = DMDestroy(&dmint);CHKERRQ(ierr); 2319 PetscFunctionReturn(0); 2320 } 2321 2322 /* step 3: construct intNodes, and intMat, and combine it with boundray data to make allNodes and allMat */ 2323 if (!tensorSpace) { 2324 if (!tensorCell) {ierr = PetscLagNodeIndicesCreateSimplexVertices(dm, &(lag->vertIndices));CHKERRQ(ierr);} 2325 2326 if (trimmed) { 2327 /* there is one dof in the interior of the a trimmed element for each full polynomial of with degree at most 2328 * order + k - dim - 1 */ 2329 if (order + PetscAbsInt(formDegree) > dim) { 2330 PetscInt sum = order + PetscAbsInt(formDegree) - dim - 1; 2331 PetscInt nDofs; 2332 2333 ierr = PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));CHKERRQ(ierr); 2334 ierr = MatGetSize(sp->intMat, &nDofs, NULL);CHKERRQ(ierr); 2335 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr); 2336 } 2337 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 2338 ierr = PetscDualSpaceCreateAllDataFromInteriorData(sp);CHKERRQ(ierr); 2339 ierr = PetscDualSpaceLagrangeCreateAllNodeIdx(sp);CHKERRQ(ierr); 2340 } else { 2341 if (!continuous) { 2342 /* if discontinuous just construct one node for each set of dofs (a set of dofs is a basis for the k-form 2343 * space) */ 2344 PetscInt sum = order; 2345 PetscInt nDofs; 2346 2347 ierr = PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));CHKERRQ(ierr); 2348 ierr = MatGetSize(sp->intMat, &nDofs, NULL);CHKERRQ(ierr); 2349 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr); 2350 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 2351 ierr = PetscObjectReference((PetscObject)(sp->intNodes));CHKERRQ(ierr); 2352 sp->allNodes = sp->intNodes; 2353 ierr = PetscObjectReference((PetscObject)(sp->intMat));CHKERRQ(ierr); 2354 sp->allMat = sp->intMat; 2355 ierr = PetscLagNodeIndicesReference(lag->intNodeIndices);CHKERRQ(ierr); 2356 lag->allNodeIndices = lag->intNodeIndices; 2357 } else { 2358 /* there is one dof in the interior of the a full element for each trimmed polynomial of with degree at most 2359 * order + k - dim, but with complementary form degree */ 2360 if (order + PetscAbsInt(formDegree) > dim) { 2361 PetscDualSpace trimmedsp; 2362 PetscDualSpace_Lag *trimmedlag; 2363 PetscQuadrature intNodes; 2364 PetscInt trFormDegree = formDegree >= 0 ? formDegree - dim : dim - PetscAbsInt(formDegree); 2365 PetscInt nDofs; 2366 Mat intMat; 2367 2368 ierr = PetscDualSpaceDuplicate(sp, &trimmedsp);CHKERRQ(ierr); 2369 ierr = PetscDualSpaceLagrangeSetTrimmed(trimmedsp, PETSC_TRUE);CHKERRQ(ierr); 2370 ierr = PetscDualSpaceSetOrder(trimmedsp, order + PetscAbsInt(formDegree) - dim);CHKERRQ(ierr); 2371 ierr = PetscDualSpaceSetFormDegree(trimmedsp, trFormDegree);CHKERRQ(ierr); 2372 trimmedlag = (PetscDualSpace_Lag *) trimmedsp->data; 2373 trimmedlag->numNodeSkip = numNodeSkip + 1; 2374 ierr = PetscDualSpaceSetUp(trimmedsp);CHKERRQ(ierr); 2375 ierr = PetscDualSpaceGetAllData(trimmedsp, &intNodes, &intMat);CHKERRQ(ierr); 2376 ierr = PetscObjectReference((PetscObject)intNodes);CHKERRQ(ierr); 2377 sp->intNodes = intNodes; 2378 ierr = PetscLagNodeIndicesReference(trimmedlag->allNodeIndices);CHKERRQ(ierr); 2379 lag->intNodeIndices = trimmedlag->allNodeIndices; 2380 ierr = PetscObjectReference((PetscObject)intMat);CHKERRQ(ierr); 2381 if (PetscAbsInt(formDegree) > 0 && PetscAbsInt(formDegree) < dim) { 2382 PetscReal *T; 2383 PetscScalar *work; 2384 PetscInt nCols, nRows; 2385 Mat intMatT; 2386 2387 ierr = MatDuplicate(intMat, MAT_COPY_VALUES, &intMatT);CHKERRQ(ierr); 2388 ierr = MatGetSize(intMat, &nRows, &nCols);CHKERRQ(ierr); 2389 ierr = PetscMalloc2(Nk * Nk, &T, nCols, &work);CHKERRQ(ierr); 2390 ierr = BiunitSimplexSymmetricFormTransformation(dim, formDegree, T);CHKERRQ(ierr); 2391 for (PetscInt row = 0; row < nRows; row++) { 2392 PetscInt nrCols; 2393 const PetscInt *rCols; 2394 const PetscScalar *rVals; 2395 2396 ierr = MatGetRow(intMat, row, &nrCols, &rCols, &rVals);CHKERRQ(ierr); 2397 if (nrCols % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in intMat matrix are not in k-form size blocks"); 2398 for (PetscInt b = 0; b < nrCols; b += Nk) { 2399 const PetscScalar *v = &rVals[b]; 2400 PetscScalar *w = &work[b]; 2401 for (PetscInt j = 0; j < Nk; j++) { 2402 w[j] = 0.; 2403 for (PetscInt i = 0; i < Nk; i++) { 2404 w[j] += v[i] * T[i * Nk + j]; 2405 } 2406 } 2407 } 2408 ierr = MatSetValuesBlocked(intMatT, 1, &row, nrCols, rCols, work, INSERT_VALUES);CHKERRQ(ierr); 2409 ierr = MatRestoreRow(intMat, row, &nrCols, &rCols, &rVals);CHKERRQ(ierr); 2410 } 2411 ierr = MatAssemblyBegin(intMatT, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2412 ierr = MatAssemblyEnd(intMatT, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2413 ierr = MatDestroy(&intMat);CHKERRQ(ierr); 2414 intMat = intMatT; 2415 ierr = PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));CHKERRQ(ierr); 2416 ierr = PetscLagNodeIndicesDuplicate(trimmedlag->allNodeIndices, &(lag->intNodeIndices));CHKERRQ(ierr); 2417 { 2418 PetscInt nNodes = lag->intNodeIndices->nNodes; 2419 PetscReal *newNodeVec = lag->intNodeIndices->nodeVec; 2420 const PetscReal *oldNodeVec = trimmedlag->allNodeIndices->nodeVec; 2421 2422 for (PetscInt n = 0; n < nNodes; n++) { 2423 PetscReal *w = &newNodeVec[n * Nk]; 2424 const PetscReal *v = &oldNodeVec[n * Nk]; 2425 2426 for (PetscInt j = 0; j < Nk; j++) { 2427 w[j] = 0.; 2428 for (PetscInt i = 0; i < Nk; i++) { 2429 w[j] += v[i] * T[i * Nk + j]; 2430 } 2431 } 2432 } 2433 } 2434 ierr = PetscFree2(T, work);CHKERRQ(ierr); 2435 } 2436 sp->intMat = intMat; 2437 ierr = MatGetSize(sp->intMat, &nDofs, NULL);CHKERRQ(ierr); 2438 ierr = PetscDualSpaceDestroy(&trimmedsp);CHKERRQ(ierr); 2439 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr); 2440 } 2441 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 2442 ierr = PetscDualSpaceCreateAllDataFromInteriorData(sp);CHKERRQ(ierr); 2443 ierr = PetscDualSpaceLagrangeCreateAllNodeIdx(sp);CHKERRQ(ierr); 2444 } 2445 } 2446 } else { 2447 PetscQuadrature intNodesTrace = NULL; 2448 PetscQuadrature intNodesFiber = NULL; 2449 PetscQuadrature intNodes = NULL; 2450 PetscLagNodeIndices intNodeIndices = NULL; 2451 Mat intMat = NULL; 2452 2453 if (PetscAbsInt(formDegree) < dim) { /* get the trace k-forms on the first facet, and the 0-forms on the edge, 2454 and wedge them together to create some of the k-form dofs */ 2455 PetscDualSpace trace, fiber; 2456 PetscDualSpace_Lag *tracel, *fiberl; 2457 Mat intMatTrace, intMatFiber; 2458 2459 if (sp->pointSpaces[tensorf]) { 2460 ierr = PetscObjectReference((PetscObject)(sp->pointSpaces[tensorf]));CHKERRQ(ierr); 2461 trace = sp->pointSpaces[tensorf]; 2462 } else { 2463 ierr = PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,formDegree,Ncopies,PETSC_TRUE,&trace);CHKERRQ(ierr); 2464 } 2465 ierr = PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,0,1,PETSC_TRUE,&fiber);CHKERRQ(ierr); 2466 tracel = (PetscDualSpace_Lag *) trace->data; 2467 fiberl = (PetscDualSpace_Lag *) fiber->data; 2468 ierr = PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));CHKERRQ(ierr); 2469 ierr = PetscDualSpaceGetInteriorData(trace, &intNodesTrace, &intMatTrace);CHKERRQ(ierr); 2470 ierr = PetscDualSpaceGetInteriorData(fiber, &intNodesFiber, &intMatFiber);CHKERRQ(ierr); 2471 if (intNodesTrace && intNodesFiber) { 2472 ierr = PetscQuadratureCreateTensor(intNodesTrace, intNodesFiber, &intNodes);CHKERRQ(ierr); 2473 ierr = MatTensorAltV(intMatTrace, intMatFiber, dim-1, formDegree, 1, 0, &intMat);CHKERRQ(ierr); 2474 ierr = PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, formDegree, fiberl->intNodeIndices, 1, 0, &intNodeIndices);CHKERRQ(ierr); 2475 } 2476 ierr = PetscObjectReference((PetscObject) intNodesTrace);CHKERRQ(ierr); 2477 ierr = PetscObjectReference((PetscObject) intNodesFiber);CHKERRQ(ierr); 2478 ierr = PetscDualSpaceDestroy(&fiber);CHKERRQ(ierr); 2479 ierr = PetscDualSpaceDestroy(&trace);CHKERRQ(ierr); 2480 } 2481 if (PetscAbsInt(formDegree) > 0) { /* get the trace (k-1)-forms on the first facet, and the 1-forms on the edge, 2482 and wedge them together to create the remaining k-form dofs */ 2483 PetscDualSpace trace, fiber; 2484 PetscDualSpace_Lag *tracel, *fiberl; 2485 PetscQuadrature intNodesTrace2, intNodesFiber2, intNodes2; 2486 PetscLagNodeIndices intNodeIndices2; 2487 Mat intMatTrace, intMatFiber, intMat2; 2488 PetscInt traceDegree = formDegree > 0 ? formDegree - 1 : formDegree + 1; 2489 PetscInt fiberDegree = formDegree > 0 ? 1 : -1; 2490 2491 ierr = PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,traceDegree,Ncopies,PETSC_TRUE,&trace);CHKERRQ(ierr); 2492 ierr = PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,fiberDegree,1,PETSC_TRUE,&fiber);CHKERRQ(ierr); 2493 tracel = (PetscDualSpace_Lag *) trace->data; 2494 fiberl = (PetscDualSpace_Lag *) fiber->data; 2495 if (!lag->vertIndices) { 2496 ierr = PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));CHKERRQ(ierr); 2497 } 2498 ierr = PetscDualSpaceGetInteriorData(trace, &intNodesTrace2, &intMatTrace);CHKERRQ(ierr); 2499 ierr = PetscDualSpaceGetInteriorData(fiber, &intNodesFiber2, &intMatFiber);CHKERRQ(ierr); 2500 if (intNodesTrace2 && intNodesFiber2) { 2501 ierr = PetscQuadratureCreateTensor(intNodesTrace2, intNodesFiber2, &intNodes2);CHKERRQ(ierr); 2502 ierr = MatTensorAltV(intMatTrace, intMatFiber, dim-1, traceDegree, 1, fiberDegree, &intMat2);CHKERRQ(ierr); 2503 ierr = PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, traceDegree, fiberl->intNodeIndices, 1, fiberDegree, &intNodeIndices2);CHKERRQ(ierr); 2504 if (!intMat) { 2505 intMat = intMat2; 2506 intNodes = intNodes2; 2507 intNodeIndices = intNodeIndices2; 2508 } else { 2509 /* merge the matrices, quadrature points, and nodes */ 2510 PetscInt nM; 2511 PetscInt nDof, nDof2; 2512 PetscInt *toMerged = NULL, *toMerged2 = NULL; 2513 PetscQuadrature merged = NULL; 2514 PetscLagNodeIndices intNodeIndicesMerged = NULL; 2515 Mat matMerged = NULL; 2516 2517 ierr = MatGetSize(intMat, &nDof, NULL);CHKERRQ(ierr); 2518 ierr = MatGetSize(intMat2, &nDof2, NULL);CHKERRQ(ierr); 2519 ierr = PetscQuadraturePointsMerge(intNodes, intNodes2, &merged, &toMerged, &toMerged2);CHKERRQ(ierr); 2520 ierr = PetscQuadratureGetData(merged, NULL, NULL, &nM, NULL, NULL);CHKERRQ(ierr); 2521 ierr = MatricesMerge(intMat, intMat2, dim, formDegree, nM, toMerged, toMerged2, &matMerged);CHKERRQ(ierr); 2522 ierr = PetscLagNodeIndicesMerge(intNodeIndices, intNodeIndices2, &intNodeIndicesMerged);CHKERRQ(ierr); 2523 ierr = PetscFree(toMerged);CHKERRQ(ierr); 2524 ierr = PetscFree(toMerged2);CHKERRQ(ierr); 2525 ierr = MatDestroy(&intMat);CHKERRQ(ierr); 2526 ierr = MatDestroy(&intMat2);CHKERRQ(ierr); 2527 ierr = PetscQuadratureDestroy(&intNodes);CHKERRQ(ierr); 2528 ierr = PetscQuadratureDestroy(&intNodes2);CHKERRQ(ierr); 2529 ierr = PetscLagNodeIndicesDestroy(&intNodeIndices);CHKERRQ(ierr); 2530 ierr = PetscLagNodeIndicesDestroy(&intNodeIndices2);CHKERRQ(ierr); 2531 intNodes = merged; 2532 intMat = matMerged; 2533 intNodeIndices = intNodeIndicesMerged; 2534 if (!trimmed) { 2535 /* I think users expect that, when a node has a full basis for the k-forms, 2536 * they should be consecutive dofs. That isn't the case for trimmed spaces, 2537 * but is for some of the nodes in untrimmed spaces, so in that case we 2538 * sort them to group them by node */ 2539 Mat intMatPerm; 2540 2541 ierr = MatPermuteByNodeIdx(intMat, intNodeIndices, &intMatPerm);CHKERRQ(ierr); 2542 ierr = MatDestroy(&intMat);CHKERRQ(ierr); 2543 intMat = intMatPerm; 2544 } 2545 } 2546 } 2547 ierr = PetscDualSpaceDestroy(&fiber);CHKERRQ(ierr); 2548 ierr = PetscDualSpaceDestroy(&trace);CHKERRQ(ierr); 2549 } 2550 ierr = PetscQuadratureDestroy(&intNodesTrace);CHKERRQ(ierr); 2551 ierr = PetscQuadratureDestroy(&intNodesFiber);CHKERRQ(ierr); 2552 sp->intNodes = intNodes; 2553 sp->intMat = intMat; 2554 lag->intNodeIndices = intNodeIndices; 2555 { 2556 PetscInt nDofs = 0; 2557 2558 if (intMat) { 2559 ierr = MatGetSize(intMat, &nDofs, NULL);CHKERRQ(ierr); 2560 } 2561 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr); 2562 } 2563 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 2564 if (continuous) { 2565 ierr = PetscDualSpaceCreateAllDataFromInteriorData(sp);CHKERRQ(ierr); 2566 ierr = PetscDualSpaceLagrangeCreateAllNodeIdx(sp);CHKERRQ(ierr); 2567 } else { 2568 ierr = PetscObjectReference((PetscObject) intNodes);CHKERRQ(ierr); 2569 sp->allNodes = intNodes; 2570 ierr = PetscObjectReference((PetscObject) intMat);CHKERRQ(ierr); 2571 sp->allMat = intMat; 2572 ierr = PetscLagNodeIndicesReference(intNodeIndices);CHKERRQ(ierr); 2573 lag->allNodeIndices = intNodeIndices; 2574 } 2575 } 2576 ierr = PetscSectionGetStorageSize(section, &sp->spdim);CHKERRQ(ierr); 2577 ierr = PetscSectionGetConstrainedStorageSize(section, &sp->spintdim);CHKERRQ(ierr); 2578 ierr = PetscDualSpaceComputeFunctionalsFromAllData(sp);CHKERRQ(ierr); 2579 ierr = PetscFree2(pStratStart, pStratEnd);CHKERRQ(ierr); 2580 ierr = DMDestroy(&dmint);CHKERRQ(ierr); 2581 PetscFunctionReturn(0); 2582 } 2583 2584 /* Create a matrix that represents the transformation that DMPlexVecGetClosure() would need 2585 * to get the representation of the dofs for a mesh point if the mesh point had this orientation 2586 * relative to the cell */ 2587 PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat) 2588 { 2589 PetscDualSpace_Lag *lag; 2590 DM dm; 2591 PetscLagNodeIndices vertIndices, intNodeIndices; 2592 PetscLagNodeIndices ni; 2593 PetscInt nodeIdxDim, nodeVecDim, nNodes; 2594 PetscInt formDegree; 2595 PetscInt *perm, *permOrnt; 2596 PetscInt *nnz; 2597 PetscInt n; 2598 PetscInt maxGroupSize; 2599 PetscScalar *V, *W, *work; 2600 Mat A; 2601 PetscErrorCode ierr; 2602 2603 PetscFunctionBegin; 2604 if (!sp->spintdim) { 2605 *symMat = NULL; 2606 PetscFunctionReturn(0); 2607 } 2608 lag = (PetscDualSpace_Lag *) sp->data; 2609 vertIndices = lag->vertIndices; 2610 intNodeIndices = lag->intNodeIndices; 2611 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 2612 ierr = PetscDualSpaceGetFormDegree(sp, &formDegree);CHKERRQ(ierr); 2613 ierr = PetscNew(&ni);CHKERRQ(ierr); 2614 ni->refct = 1; 2615 ni->nodeIdxDim = nodeIdxDim = intNodeIndices->nodeIdxDim; 2616 ni->nodeVecDim = nodeVecDim = intNodeIndices->nodeVecDim; 2617 ni->nNodes = nNodes = intNodeIndices->nNodes; 2618 ierr = PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));CHKERRQ(ierr); 2619 ierr = PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));CHKERRQ(ierr); 2620 /* push forward the dofs by the symmetry of the reference element induced by ornt */ 2621 ierr = PetscLagNodeIndicesPushForward(dm, vertIndices, 0, vertIndices, intNodeIndices, ornt, formDegree, ni->nodeIdx, ni->nodeVec);CHKERRQ(ierr); 2622 /* get the revlex order for both the original and transformed dofs */ 2623 ierr = PetscLagNodeIndicesGetPermutation(intNodeIndices, &perm);CHKERRQ(ierr); 2624 ierr = PetscLagNodeIndicesGetPermutation(ni, &permOrnt);CHKERRQ(ierr); 2625 ierr = PetscMalloc1(nNodes, &nnz);CHKERRQ(ierr); 2626 for (n = 0, maxGroupSize = 0; n < nNodes;) { /* incremented in the loop */ 2627 PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]); 2628 PetscInt m, nEnd; 2629 PetscInt groupSize; 2630 /* for each group of dofs that have the same nodeIdx coordinate */ 2631 for (nEnd = n + 1; nEnd < nNodes; nEnd++) { 2632 PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]); 2633 PetscInt d; 2634 2635 /* compare the oriented permutation indices */ 2636 for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break; 2637 if (d < nodeIdxDim) break; 2638 } 2639 /* permOrnt[[n, nEnd)] is a group of dofs that, under the symmetry are at the same location */ 2640 2641 /* the symmetry had better map the group of dofs with the same permuted nodeIdx 2642 * to a group of dofs with the same size, otherwise we messed up */ 2643 if (PetscDefined(USE_DEBUG)) { 2644 PetscInt m; 2645 PetscInt *nind = &(intNodeIndices->nodeIdx[perm[n] * nodeIdxDim]); 2646 2647 for (m = n + 1; m < nEnd; m++) { 2648 PetscInt *mind = &(intNodeIndices->nodeIdx[perm[m] * nodeIdxDim]); 2649 PetscInt d; 2650 2651 /* compare the oriented permutation indices */ 2652 for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break; 2653 if (d < nodeIdxDim) break; 2654 } 2655 if (m < nEnd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs with same index after symmetry not same block size"); 2656 } 2657 groupSize = nEnd - n; 2658 /* each pushforward dof vector will be expressed in a basis of the unpermuted dofs */ 2659 for (m = n; m < nEnd; m++) nnz[permOrnt[m]] = groupSize; 2660 2661 maxGroupSize = PetscMax(maxGroupSize, nEnd - n); 2662 n = nEnd; 2663 } 2664 if (maxGroupSize > nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs are not in blocks that can be solved"); 2665 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes, nNodes, 0, nnz, &A);CHKERRQ(ierr); 2666 ierr = PetscFree(nnz);CHKERRQ(ierr); 2667 ierr = PetscMalloc3(maxGroupSize * nodeVecDim, &V, maxGroupSize * nodeVecDim, &W, nodeVecDim * 2, &work);CHKERRQ(ierr); 2668 for (n = 0; n < nNodes;) { /* incremented in the loop */ 2669 PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]); 2670 PetscInt nEnd; 2671 PetscInt m; 2672 PetscInt groupSize; 2673 for (nEnd = n + 1; nEnd < nNodes; nEnd++) { 2674 PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]); 2675 PetscInt d; 2676 2677 /* compare the oriented permutation indices */ 2678 for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break; 2679 if (d < nodeIdxDim) break; 2680 } 2681 groupSize = nEnd - n; 2682 /* get all of the vectors from the original and all of the pushforward vectors */ 2683 for (m = n; m < nEnd; m++) { 2684 PetscInt d; 2685 2686 for (d = 0; d < nodeVecDim; d++) { 2687 V[(m - n) * nodeVecDim + d] = intNodeIndices->nodeVec[perm[m] * nodeVecDim + d]; 2688 W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d]; 2689 } 2690 } 2691 /* now we have to solve for W in terms of V: the systems isn't always square, but the span 2692 * of V and W should always be the same, so the solution of the normal equations works */ 2693 { 2694 char transpose = 'N'; 2695 PetscBLASInt bm = nodeVecDim; 2696 PetscBLASInt bn = groupSize; 2697 PetscBLASInt bnrhs = groupSize; 2698 PetscBLASInt blda = bm; 2699 PetscBLASInt bldb = bm; 2700 PetscBLASInt blwork = 2 * nodeVecDim; 2701 PetscBLASInt info; 2702 2703 PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&bm,&bn,&bnrhs,V,&blda,W,&bldb,work,&blwork, &info)); 2704 if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS"); 2705 /* repack */ 2706 { 2707 PetscInt i, j; 2708 2709 for (i = 0; i < groupSize; i++) { 2710 for (j = 0; j < groupSize; j++) { 2711 /* notice the different leading dimension */ 2712 V[i * groupSize + j] = W[i * nodeVecDim + j]; 2713 } 2714 } 2715 } 2716 if (PetscDefined(USE_DEBUG)) { 2717 PetscReal res; 2718 2719 /* check that the normal error is 0 */ 2720 for (m = n; m < nEnd; m++) { 2721 PetscInt d; 2722 2723 for (d = 0; d < nodeVecDim; d++) { 2724 W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d]; 2725 } 2726 } 2727 res = 0.; 2728 for (PetscInt i = 0; i < groupSize; i++) { 2729 for (PetscInt j = 0; j < nodeVecDim; j++) { 2730 for (PetscInt k = 0; k < groupSize; k++) { 2731 W[i * nodeVecDim + j] -= V[i * groupSize + k] * intNodeIndices->nodeVec[perm[n+k] * nodeVecDim + j]; 2732 } 2733 res += PetscAbsScalar(W[i * nodeVecDim + j]); 2734 } 2735 } 2736 if (res > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Dof block did not solve"); 2737 } 2738 } 2739 ierr = MatSetValues(A, groupSize, &permOrnt[n], groupSize, &perm[n], V, INSERT_VALUES);CHKERRQ(ierr); 2740 n = nEnd; 2741 } 2742 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2743 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2744 *symMat = A; 2745 ierr = PetscFree3(V,W,work);CHKERRQ(ierr); 2746 ierr = PetscLagNodeIndicesDestroy(&ni);CHKERRQ(ierr); 2747 PetscFunctionReturn(0); 2748 } 2749 2750 #define BaryIndex(perEdge,a,b,c) (((b)*(2*perEdge+1-(b)))/2)+(c) 2751 2752 #define CartIndex(perEdge,a,b) (perEdge*(a)+b) 2753 2754 /* the existing interface for symmetries is insufficient for all cases: 2755 * - it should be sufficient for form degrees that are scalar (0 and n) 2756 * - it should be sufficient for hypercube dofs 2757 * - it isn't sufficient for simplex cells with non-scalar form degrees if 2758 * there are any dofs in the interior 2759 * 2760 * We compute the general transformation matrices, and if they fit, we return them, 2761 * otherwise we error (but we should probably change the interface to allow for 2762 * these symmetries) 2763 */ 2764 static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips) 2765 { 2766 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 2767 PetscInt dim, order, Nc; 2768 PetscErrorCode ierr; 2769 2770 PetscFunctionBegin; 2771 ierr = PetscDualSpaceGetOrder(sp,&order);CHKERRQ(ierr); 2772 ierr = PetscDualSpaceGetNumComponents(sp,&Nc);CHKERRQ(ierr); 2773 ierr = DMGetDimension(sp->dm,&dim);CHKERRQ(ierr); 2774 if (!lag->symComputed) { /* store symmetries */ 2775 PetscInt pStart, pEnd, p; 2776 PetscInt numPoints; 2777 PetscInt numFaces; 2778 PetscInt spintdim; 2779 PetscInt ***symperms; 2780 PetscScalar ***symflips; 2781 2782 ierr = DMPlexGetChart(sp->dm, &pStart, &pEnd);CHKERRQ(ierr); 2783 numPoints = pEnd - pStart; 2784 { 2785 DMPolytopeType ct; 2786 /* The number of arrangements is no longer based on the number of faces */ 2787 ierr = DMPlexGetCellType(sp->dm, 0, &ct);CHKERRQ(ierr); 2788 numFaces = DMPolytopeTypeGetNumArrangments(ct) / 2; 2789 } 2790 ierr = PetscCalloc1(numPoints,&symperms);CHKERRQ(ierr); 2791 ierr = PetscCalloc1(numPoints,&symflips);CHKERRQ(ierr); 2792 spintdim = sp->spintdim; 2793 /* The nodal symmetry behavior is not present when tensorSpace != tensorCell: someone might want this for the "S" 2794 * family of FEEC spaces. Most used in particular are discontinuous polynomial L2 spaces in tensor cells, where 2795 * the symmetries are not necessary for FE assembly. So for now we assume this is the case and don't return 2796 * symmetries if tensorSpace != tensorCell */ 2797 if (spintdim && 0 < dim && dim < 3 && (lag->tensorSpace == lag->tensorCell)) { /* compute self symmetries */ 2798 PetscInt **cellSymperms; 2799 PetscScalar **cellSymflips; 2800 PetscInt ornt; 2801 PetscInt nCopies = Nc / lag->intNodeIndices->nodeVecDim; 2802 PetscInt nNodes = lag->intNodeIndices->nNodes; 2803 2804 lag->numSelfSym = 2 * numFaces; 2805 lag->selfSymOff = numFaces; 2806 ierr = PetscCalloc1(2*numFaces,&cellSymperms);CHKERRQ(ierr); 2807 ierr = PetscCalloc1(2*numFaces,&cellSymflips);CHKERRQ(ierr); 2808 /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */ 2809 symperms[0] = &cellSymperms[numFaces]; 2810 symflips[0] = &cellSymflips[numFaces]; 2811 if (lag->intNodeIndices->nodeVecDim * nCopies != Nc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs"); 2812 if (nNodes * nCopies != spintdim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs"); 2813 for (ornt = -numFaces; ornt < numFaces; ornt++) { /* for every symmetry, compute the symmetry matrix, and extract rows to see if it fits in the perm + flip framework */ 2814 Mat symMat; 2815 PetscInt *perm; 2816 PetscScalar *flips; 2817 PetscInt i; 2818 2819 if (!ornt) continue; 2820 ierr = PetscMalloc1(spintdim, &perm);CHKERRQ(ierr); 2821 ierr = PetscCalloc1(spintdim, &flips);CHKERRQ(ierr); 2822 for (i = 0; i < spintdim; i++) perm[i] = -1; 2823 ierr = PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, ornt, &symMat);CHKERRQ(ierr); 2824 for (i = 0; i < nNodes; i++) { 2825 PetscInt ncols; 2826 PetscInt j, k; 2827 const PetscInt *cols; 2828 const PetscScalar *vals; 2829 PetscBool nz_seen = PETSC_FALSE; 2830 2831 ierr = MatGetRow(symMat, i, &ncols, &cols, &vals);CHKERRQ(ierr); 2832 for (j = 0; j < ncols; j++) { 2833 if (PetscAbsScalar(vals[j]) > PETSC_SMALL) { 2834 if (nz_seen) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips"); 2835 nz_seen = PETSC_TRUE; 2836 if (PetscAbsReal(PetscAbsScalar(vals[j]) - PetscRealConstant(1.)) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips"); 2837 if (PetscAbsReal(PetscImaginaryPart(vals[j])) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips"); 2838 if (perm[cols[j] * nCopies] >= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips"); 2839 for (k = 0; k < nCopies; k++) { 2840 perm[cols[j] * nCopies + k] = i * nCopies + k; 2841 } 2842 if (PetscRealPart(vals[j]) < 0.) { 2843 for (k = 0; k < nCopies; k++) { 2844 flips[i * nCopies + k] = -1.; 2845 } 2846 } else { 2847 for (k = 0; k < nCopies; k++) { 2848 flips[i * nCopies + k] = 1.; 2849 } 2850 } 2851 } 2852 } 2853 ierr = MatRestoreRow(symMat, i, &ncols, &cols, &vals);CHKERRQ(ierr); 2854 } 2855 ierr = MatDestroy(&symMat);CHKERRQ(ierr); 2856 /* if there were no sign flips, keep NULL */ 2857 for (i = 0; i < spintdim; i++) if (flips[i] != 1.) break; 2858 if (i == spintdim) { 2859 ierr = PetscFree(flips);CHKERRQ(ierr); 2860 flips = NULL; 2861 } 2862 /* if the permutation is identity, keep NULL */ 2863 for (i = 0; i < spintdim; i++) if (perm[i] != i) break; 2864 if (i == spintdim) { 2865 ierr = PetscFree(perm);CHKERRQ(ierr); 2866 perm = NULL; 2867 } 2868 symperms[0][ornt] = perm; 2869 symflips[0][ornt] = flips; 2870 } 2871 /* if no orientations produced non-identity permutations, keep NULL */ 2872 for (ornt = -numFaces; ornt < numFaces; ornt++) if (symperms[0][ornt]) break; 2873 if (ornt == numFaces) { 2874 ierr = PetscFree(cellSymperms);CHKERRQ(ierr); 2875 symperms[0] = NULL; 2876 } 2877 /* if no orientations produced sign flips, keep NULL */ 2878 for (ornt = -numFaces; ornt < numFaces; ornt++) if (symflips[0][ornt]) break; 2879 if (ornt == numFaces) { 2880 ierr = PetscFree(cellSymflips);CHKERRQ(ierr); 2881 symflips[0] = NULL; 2882 } 2883 } 2884 { /* get the symmetries of closure points */ 2885 PetscInt closureSize = 0; 2886 PetscInt *closure = NULL; 2887 PetscInt r; 2888 2889 ierr = DMPlexGetTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2890 for (r = 0; r < closureSize; r++) { 2891 PetscDualSpace psp; 2892 PetscInt point = closure[2 * r]; 2893 PetscInt pspintdim; 2894 const PetscInt ***psymperms = NULL; 2895 const PetscScalar ***psymflips = NULL; 2896 2897 if (!point) continue; 2898 ierr = PetscDualSpaceGetPointSubspace(sp, point, &psp);CHKERRQ(ierr); 2899 if (!psp) continue; 2900 ierr = PetscDualSpaceGetInteriorDimension(psp, &pspintdim);CHKERRQ(ierr); 2901 if (!pspintdim) continue; 2902 ierr = PetscDualSpaceGetSymmetries(psp,&psymperms,&psymflips);CHKERRQ(ierr); 2903 symperms[r] = (PetscInt **) (psymperms ? psymperms[0] : NULL); 2904 symflips[r] = (PetscScalar **) (psymflips ? psymflips[0] : NULL); 2905 } 2906 ierr = DMPlexRestoreTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2907 } 2908 for (p = 0; p < pEnd; p++) if (symperms[p]) break; 2909 if (p == pEnd) { 2910 ierr = PetscFree(symperms);CHKERRQ(ierr); 2911 symperms = NULL; 2912 } 2913 for (p = 0; p < pEnd; p++) if (symflips[p]) break; 2914 if (p == pEnd) { 2915 ierr = PetscFree(symflips);CHKERRQ(ierr); 2916 symflips = NULL; 2917 } 2918 lag->symperms = symperms; 2919 lag->symflips = symflips; 2920 lag->symComputed = PETSC_TRUE; 2921 } 2922 if (perms) *perms = (const PetscInt ***) lag->symperms; 2923 if (flips) *flips = (const PetscScalar ***) lag->symflips; 2924 PetscFunctionReturn(0); 2925 } 2926 2927 static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous) 2928 { 2929 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 2930 2931 PetscFunctionBegin; 2932 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2933 PetscValidPointer(continuous, 2); 2934 *continuous = lag->continuous; 2935 PetscFunctionReturn(0); 2936 } 2937 2938 static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous) 2939 { 2940 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 2941 2942 PetscFunctionBegin; 2943 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2944 lag->continuous = continuous; 2945 PetscFunctionReturn(0); 2946 } 2947 2948 /*@ 2949 PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity 2950 2951 Not Collective 2952 2953 Input Parameter: 2954 . sp - the PetscDualSpace 2955 2956 Output Parameter: 2957 . continuous - flag for element continuity 2958 2959 Level: intermediate 2960 2961 .seealso: PetscDualSpaceLagrangeSetContinuity() 2962 @*/ 2963 PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous) 2964 { 2965 PetscErrorCode ierr; 2966 2967 PetscFunctionBegin; 2968 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2969 PetscValidPointer(continuous, 2); 2970 ierr = PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));CHKERRQ(ierr); 2971 PetscFunctionReturn(0); 2972 } 2973 2974 /*@ 2975 PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous 2976 2977 Logically Collective on sp 2978 2979 Input Parameters: 2980 + sp - the PetscDualSpace 2981 - continuous - flag for element continuity 2982 2983 Options Database: 2984 . -petscdualspace_lagrange_continuity <bool> 2985 2986 Level: intermediate 2987 2988 .seealso: PetscDualSpaceLagrangeGetContinuity() 2989 @*/ 2990 PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous) 2991 { 2992 PetscErrorCode ierr; 2993 2994 PetscFunctionBegin; 2995 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2996 PetscValidLogicalCollectiveBool(sp, continuous, 2); 2997 ierr = PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));CHKERRQ(ierr); 2998 PetscFunctionReturn(0); 2999 } 3000 3001 static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor) 3002 { 3003 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 3004 3005 PetscFunctionBegin; 3006 *tensor = lag->tensorSpace; 3007 PetscFunctionReturn(0); 3008 } 3009 3010 static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor) 3011 { 3012 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 3013 3014 PetscFunctionBegin; 3015 lag->tensorSpace = tensor; 3016 PetscFunctionReturn(0); 3017 } 3018 3019 static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Lagrange(PetscDualSpace sp, PetscBool *trimmed) 3020 { 3021 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 3022 3023 PetscFunctionBegin; 3024 *trimmed = lag->trimmed; 3025 PetscFunctionReturn(0); 3026 } 3027 3028 static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp, PetscBool trimmed) 3029 { 3030 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 3031 3032 PetscFunctionBegin; 3033 lag->trimmed = trimmed; 3034 PetscFunctionReturn(0); 3035 } 3036 3037 static PetscErrorCode PetscDualSpaceLagrangeGetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent) 3038 { 3039 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 3040 3041 PetscFunctionBegin; 3042 if (nodeType) *nodeType = lag->nodeType; 3043 if (boundary) *boundary = lag->endNodes; 3044 if (exponent) *exponent = lag->nodeExponent; 3045 PetscFunctionReturn(0); 3046 } 3047 3048 static PetscErrorCode PetscDualSpaceLagrangeSetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent) 3049 { 3050 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 3051 3052 PetscFunctionBegin; 3053 if (nodeType == PETSCDTNODES_GAUSSJACOBI && exponent <= -1.) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Exponent must be > -1"); 3054 lag->nodeType = nodeType; 3055 lag->endNodes = boundary; 3056 lag->nodeExponent = exponent; 3057 PetscFunctionReturn(0); 3058 } 3059 3060 static PetscErrorCode PetscDualSpaceLagrangeGetUseMoments_Lagrange(PetscDualSpace sp, PetscBool *useMoments) 3061 { 3062 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 3063 3064 PetscFunctionBegin; 3065 *useMoments = lag->useMoments; 3066 PetscFunctionReturn(0); 3067 } 3068 3069 static PetscErrorCode PetscDualSpaceLagrangeSetUseMoments_Lagrange(PetscDualSpace sp, PetscBool useMoments) 3070 { 3071 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 3072 3073 PetscFunctionBegin; 3074 lag->useMoments = useMoments; 3075 PetscFunctionReturn(0); 3076 } 3077 3078 static PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt *momentOrder) 3079 { 3080 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 3081 3082 PetscFunctionBegin; 3083 *momentOrder = lag->momentOrder; 3084 PetscFunctionReturn(0); 3085 } 3086 3087 static PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt momentOrder) 3088 { 3089 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 3090 3091 PetscFunctionBegin; 3092 lag->momentOrder = momentOrder; 3093 PetscFunctionReturn(0); 3094 } 3095 3096 /*@ 3097 PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space 3098 3099 Not collective 3100 3101 Input Parameter: 3102 . sp - The PetscDualSpace 3103 3104 Output Parameter: 3105 . tensor - Whether the dual space has tensor layout (vs. simplicial) 3106 3107 Level: intermediate 3108 3109 .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate() 3110 @*/ 3111 PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor) 3112 { 3113 PetscErrorCode ierr; 3114 3115 PetscFunctionBegin; 3116 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3117 PetscValidPointer(tensor, 2); 3118 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));CHKERRQ(ierr); 3119 PetscFunctionReturn(0); 3120 } 3121 3122 /*@ 3123 PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space 3124 3125 Not collective 3126 3127 Input Parameters: 3128 + sp - The PetscDualSpace 3129 - tensor - Whether the dual space has tensor layout (vs. simplicial) 3130 3131 Level: intermediate 3132 3133 .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate() 3134 @*/ 3135 PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor) 3136 { 3137 PetscErrorCode ierr; 3138 3139 PetscFunctionBegin; 3140 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3141 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));CHKERRQ(ierr); 3142 PetscFunctionReturn(0); 3143 } 3144 3145 /*@ 3146 PetscDualSpaceLagrangeGetTrimmed - Get the trimmed nature of the dual space 3147 3148 Not collective 3149 3150 Input Parameter: 3151 . sp - The PetscDualSpace 3152 3153 Output Parameter: 3154 . trimmed - Whether the dual space represents to dual basis of a trimmed polynomial space (e.g. Raviart-Thomas and higher order / other form degree variants) 3155 3156 Level: intermediate 3157 3158 .seealso: PetscDualSpaceLagrangeSetTrimmed(), PetscDualSpaceCreate() 3159 @*/ 3160 PetscErrorCode PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp, PetscBool *trimmed) 3161 { 3162 PetscErrorCode ierr; 3163 3164 PetscFunctionBegin; 3165 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3166 PetscValidPointer(trimmed, 2); 3167 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTrimmed_C",(PetscDualSpace,PetscBool *),(sp,trimmed));CHKERRQ(ierr); 3168 PetscFunctionReturn(0); 3169 } 3170 3171 /*@ 3172 PetscDualSpaceLagrangeSetTrimmed - Set the trimmed nature of the dual space 3173 3174 Not collective 3175 3176 Input Parameters: 3177 + sp - The PetscDualSpace 3178 - trimmed - Whether the dual space represents to dual basis of a trimmed polynomial space (e.g. Raviart-Thomas and higher order / other form degree variants) 3179 3180 Level: intermediate 3181 3182 .seealso: PetscDualSpaceLagrangeGetTrimmed(), PetscDualSpaceCreate() 3183 @*/ 3184 PetscErrorCode PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp, PetscBool trimmed) 3185 { 3186 PetscErrorCode ierr; 3187 3188 PetscFunctionBegin; 3189 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3190 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTrimmed_C",(PetscDualSpace,PetscBool),(sp,trimmed));CHKERRQ(ierr); 3191 PetscFunctionReturn(0); 3192 } 3193 3194 /*@ 3195 PetscDualSpaceLagrangeGetNodeType - Get a description of how nodes are laid out for Lagrange polynomials in this 3196 dual space 3197 3198 Not collective 3199 3200 Input Parameter: 3201 . sp - The PetscDualSpace 3202 3203 Output Parameters: 3204 + nodeType - The type of nodes 3205 . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that 3206 include the boundary are Gauss-Lobatto-Jacobi nodes) 3207 - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function 3208 '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type 3209 3210 Level: advanced 3211 3212 .seealso: PetscDTNodeType, PetscDualSpaceLagrangeSetNodeType() 3213 @*/ 3214 PetscErrorCode PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent) 3215 { 3216 PetscErrorCode ierr; 3217 3218 PetscFunctionBegin; 3219 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3220 if (nodeType) PetscValidPointer(nodeType, 2); 3221 if (boundary) PetscValidPointer(boundary, 3); 3222 if (exponent) PetscValidPointer(exponent, 4); 3223 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeGetNodeType_C",(PetscDualSpace,PetscDTNodeType *,PetscBool *,PetscReal *),(sp,nodeType,boundary,exponent));CHKERRQ(ierr); 3224 PetscFunctionReturn(0); 3225 } 3226 3227 /*@ 3228 PetscDualSpaceLagrangeSetNodeType - Set a description of how nodes are laid out for Lagrange polynomials in this 3229 dual space 3230 3231 Logically collective 3232 3233 Input Parameters: 3234 + sp - The PetscDualSpace 3235 . nodeType - The type of nodes 3236 . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that 3237 include the boundary are Gauss-Lobatto-Jacobi nodes) 3238 - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function 3239 '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type 3240 3241 Level: advanced 3242 3243 .seealso: PetscDTNodeType, PetscDualSpaceLagrangeGetNodeType() 3244 @*/ 3245 PetscErrorCode PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent) 3246 { 3247 PetscErrorCode ierr; 3248 3249 PetscFunctionBegin; 3250 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3251 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetNodeType_C",(PetscDualSpace,PetscDTNodeType,PetscBool,PetscReal),(sp,nodeType,boundary,exponent));CHKERRQ(ierr); 3252 PetscFunctionReturn(0); 3253 } 3254 3255 /*@ 3256 PetscDualSpaceLagrangeGetUseMoments - Get the flag for using moment functionals 3257 3258 Not collective 3259 3260 Input Parameter: 3261 . sp - The PetscDualSpace 3262 3263 Output Parameter: 3264 . useMoments - Moment flag 3265 3266 Level: advanced 3267 3268 .seealso: PetscDualSpaceLagrangeSetUseMoments() 3269 @*/ 3270 PetscErrorCode PetscDualSpaceLagrangeGetUseMoments(PetscDualSpace sp, PetscBool *useMoments) 3271 { 3272 PetscErrorCode ierr; 3273 3274 PetscFunctionBegin; 3275 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3276 PetscValidBoolPointer(useMoments, 2); 3277 ierr = PetscUseMethod(sp,"PetscDualSpaceLagrangeGetUseMoments_C",(PetscDualSpace,PetscBool *),(sp,useMoments));CHKERRQ(ierr); 3278 PetscFunctionReturn(0); 3279 } 3280 3281 /*@ 3282 PetscDualSpaceLagrangeSetUseMoments - Set the flag for moment functionals 3283 3284 Logically collective 3285 3286 Input Parameters: 3287 + sp - The PetscDualSpace 3288 - useMoments - The flag for moment functionals 3289 3290 Level: advanced 3291 3292 .seealso: PetscDualSpaceLagrangeGetUseMoments() 3293 @*/ 3294 PetscErrorCode PetscDualSpaceLagrangeSetUseMoments(PetscDualSpace sp, PetscBool useMoments) 3295 { 3296 PetscErrorCode ierr; 3297 3298 PetscFunctionBegin; 3299 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3300 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetUseMoments_C",(PetscDualSpace,PetscBool),(sp,useMoments));CHKERRQ(ierr); 3301 PetscFunctionReturn(0); 3302 } 3303 3304 /*@ 3305 PetscDualSpaceLagrangeGetMomentOrder - Get the order for moment integration 3306 3307 Not collective 3308 3309 Input Parameter: 3310 . sp - The PetscDualSpace 3311 3312 Output Parameter: 3313 . order - Moment integration order 3314 3315 Level: advanced 3316 3317 .seealso: PetscDualSpaceLagrangeSetMomentOrder() 3318 @*/ 3319 PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder(PetscDualSpace sp, PetscInt *order) 3320 { 3321 PetscErrorCode ierr; 3322 3323 PetscFunctionBegin; 3324 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3325 PetscValidIntPointer(order, 2); 3326 ierr = PetscUseMethod(sp,"PetscDualSpaceLagrangeGetMomentOrder_C",(PetscDualSpace,PetscInt *),(sp,order));CHKERRQ(ierr); 3327 PetscFunctionReturn(0); 3328 } 3329 3330 /*@ 3331 PetscDualSpaceLagrangeSetMomentOrder - Set the order for moment integration 3332 3333 Logically collective 3334 3335 Input Parameters: 3336 + sp - The PetscDualSpace 3337 - order - The order for moment integration 3338 3339 Level: advanced 3340 3341 .seealso: PetscDualSpaceLagrangeGetMomentOrder() 3342 @*/ 3343 PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder(PetscDualSpace sp, PetscInt order) 3344 { 3345 PetscErrorCode ierr; 3346 3347 PetscFunctionBegin; 3348 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3349 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetMomentOrder_C",(PetscDualSpace,PetscInt),(sp,order));CHKERRQ(ierr); 3350 PetscFunctionReturn(0); 3351 } 3352 3353 static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp) 3354 { 3355 PetscFunctionBegin; 3356 sp->ops->destroy = PetscDualSpaceDestroy_Lagrange; 3357 sp->ops->view = PetscDualSpaceView_Lagrange; 3358 sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange; 3359 sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange; 3360 sp->ops->setup = PetscDualSpaceSetUp_Lagrange; 3361 sp->ops->createheightsubspace = NULL; 3362 sp->ops->createpointsubspace = NULL; 3363 sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Lagrange; 3364 sp->ops->apply = PetscDualSpaceApplyDefault; 3365 sp->ops->applyall = PetscDualSpaceApplyAllDefault; 3366 sp->ops->applyint = PetscDualSpaceApplyInteriorDefault; 3367 sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault; 3368 sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault; 3369 PetscFunctionReturn(0); 3370 } 3371 3372 /*MC 3373 PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals 3374 3375 Level: intermediate 3376 3377 .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType() 3378 M*/ 3379 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp) 3380 { 3381 PetscDualSpace_Lag *lag; 3382 PetscErrorCode ierr; 3383 3384 PetscFunctionBegin; 3385 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 3386 ierr = PetscNewLog(sp,&lag);CHKERRQ(ierr); 3387 sp->data = lag; 3388 3389 lag->tensorCell = PETSC_FALSE; 3390 lag->tensorSpace = PETSC_FALSE; 3391 lag->continuous = PETSC_TRUE; 3392 lag->numCopies = PETSC_DEFAULT; 3393 lag->numNodeSkip = PETSC_DEFAULT; 3394 lag->nodeType = PETSCDTNODES_DEFAULT; 3395 lag->useMoments = PETSC_FALSE; 3396 lag->momentOrder = 0; 3397 3398 ierr = PetscDualSpaceInitialize_Lagrange(sp);CHKERRQ(ierr); 3399 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);CHKERRQ(ierr); 3400 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);CHKERRQ(ierr); 3401 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);CHKERRQ(ierr); 3402 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);CHKERRQ(ierr); 3403 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Lagrange);CHKERRQ(ierr); 3404 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Lagrange);CHKERRQ(ierr); 3405 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Lagrange);CHKERRQ(ierr); 3406 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Lagrange);CHKERRQ(ierr); 3407 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetUseMoments_C", PetscDualSpaceLagrangeGetUseMoments_Lagrange);CHKERRQ(ierr); 3408 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetUseMoments_C", PetscDualSpaceLagrangeSetUseMoments_Lagrange);CHKERRQ(ierr); 3409 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetMomentOrder_C", PetscDualSpaceLagrangeGetMomentOrder_Lagrange);CHKERRQ(ierr); 3410 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetMomentOrder_C", PetscDualSpaceLagrangeSetMomentOrder_Lagrange);CHKERRQ(ierr); 3411 PetscFunctionReturn(0); 3412 } 3413