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