1 #include <petsc/private/dmfieldimpl.h> /*I "petscdmfield.h" I*/ 2 #include <petscfe.h> 3 #include <petscdmplex.h> 4 #include <petscds.h> 5 6 typedef struct _n_DMField_DS 7 { 8 PetscInt fieldNum; 9 Vec vec; 10 PetscInt height; 11 PetscObject *disc; 12 PetscBool multifieldVec; 13 } 14 DMField_DS; 15 16 static PetscErrorCode DMFieldDestroy_DS(DMField field) 17 { 18 DMField_DS *dsfield; 19 PetscInt i; 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 dsfield = (DMField_DS *) field->data; 24 ierr = VecDestroy(&dsfield->vec);CHKERRQ(ierr); 25 for (i = 0; i < dsfield->height; i++) { 26 ierr = PetscObjectDereference(dsfield->disc[i]);CHKERRQ(ierr); 27 } 28 ierr = PetscFree(dsfield->disc);CHKERRQ(ierr); 29 ierr = PetscFree(dsfield);CHKERRQ(ierr); 30 PetscFunctionReturn(0); 31 } 32 33 static PetscErrorCode DMFieldView_DS(DMField field,PetscViewer viewer) 34 { 35 DMField_DS *dsfield = (DMField_DS *) field->data; 36 PetscBool iascii; 37 PetscObject disc; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 42 disc = dsfield->disc[0]; 43 if (iascii) { 44 PetscViewerASCIIPrintf(viewer, "PetscDS field %D\n",dsfield->fieldNum);CHKERRQ(ierr); 45 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 46 ierr = PetscObjectView(disc,viewer);CHKERRQ(ierr); 47 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 48 } 49 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 50 if (dsfield->multifieldVec) { 51 SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"View of subfield not implemented yet"); 52 } else { 53 ierr = VecView(dsfield->vec,viewer);CHKERRQ(ierr); 54 } 55 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 56 PetscFunctionReturn(0); 57 } 58 59 static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H) 60 { 61 PetscFunctionBegin; 62 SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented yet"); 63 PetscFunctionReturn(0); 64 } 65 66 static PetscErrorCode DMFieldDSGetHeightDisc(DMField field, PetscInt height, PetscObject *disc) 67 { 68 DMField_DS *dsfield = (DMField_DS *) field->data; 69 PetscErrorCode ierr; 70 71 PetscFunctionBegin; 72 if (!dsfield->disc[height]) { 73 PetscClassId id; 74 75 ierr = PetscObjectGetClassId(dsfield->disc[0],&id);CHKERRQ(ierr); 76 if (id == PETSCFE_CLASSID) { 77 PetscFE fe = (PetscFE) dsfield->disc[0]; 78 79 ierr = PetscFECreateHeightTrace(fe,height,(PetscFE *)&dsfield->disc[height]);CHKERRQ(ierr); 80 } 81 } 82 *disc = dsfield->disc[height]; 83 PetscFunctionReturn(0); 84 } 85 86 #define DMFieldDSdot(y,A,b,m,n,c,cast) \ 87 do { \ 88 PetscInt _i, _j, _k; \ 89 for (_i = 0; _i < (m); _i++) { \ 90 for (_k = 0; _k < (c); _k++) { \ 91 (y)[_i * (c) + _k] = 0.; \ 92 } \ 93 for (_j = 0; _j < (n); _j++) { \ 94 for (_k = 0; _k < (c); _k++) { \ 95 (y)[_i * (c) + _k] += (A)[(_i * (n) + _j) * (c) + _k] * cast((b)[_j]); \ 96 } \ 97 } \ 98 } \ 99 } while (0) 100 101 static PetscErrorCode DMFieldEvaluateFE_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H) 102 { 103 DMField_DS *dsfield = (DMField_DS *) field->data; 104 DM dm; 105 PetscObject disc; 106 PetscClassId classid; 107 PetscInt nq, nc, dim, meshDim, numCells; 108 PetscSection section; 109 const PetscReal *qpoints; 110 PetscBool isStride; 111 const PetscInt *points = NULL; 112 PetscInt sfirst = -1, stride = -1; 113 PetscErrorCode ierr; 114 115 PetscFunctionBeginHot; 116 dm = field->dm; 117 nc = field->numComponents; 118 ierr = PetscQuadratureGetData(quad,&dim,NULL,&nq,&qpoints,NULL);CHKERRQ(ierr); 119 ierr = DMFieldDSGetHeightDisc(field,dsfield->height - 1 - dim,&disc);CHKERRQ(ierr); 120 ierr = DMGetDimension(dm,&meshDim);CHKERRQ(ierr); 121 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 122 ierr = PetscSectionGetField(section,dsfield->fieldNum,§ion);CHKERRQ(ierr); 123 ierr = PetscObjectGetClassId(disc,&classid);CHKERRQ(ierr); 124 /* TODO: batch */ 125 ierr = PetscObjectTypeCompare((PetscObject)pointIS,ISSTRIDE,&isStride);CHKERRQ(ierr); 126 ierr = ISGetLocalSize(pointIS,&numCells);CHKERRQ(ierr); 127 if (isStride) { 128 ierr = ISStrideGetInfo(pointIS,&sfirst,&stride);CHKERRQ(ierr); 129 } else { 130 ierr = ISGetIndices(pointIS,&points);CHKERRQ(ierr); 131 } 132 if (classid == PETSCFE_CLASSID) { 133 PetscFE fe = (PetscFE) disc; 134 PetscInt feDim, i; 135 PetscReal *fB = NULL, *fD = NULL, *fH = NULL; 136 137 if (dim == meshDim - 1) { 138 /* TODO */ 139 } 140 ierr = PetscFEGetDimension(fe,&feDim);CHKERRQ(ierr); 141 ierr = PetscFEGetTabulation(fe,nq,qpoints,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);CHKERRQ(ierr); 142 for (i = 0; i < numCells; i++) { 143 PetscInt c = isStride ? (sfirst + i * stride) : points[i]; 144 PetscInt closureSize; 145 PetscScalar *elem = NULL; 146 147 ierr = DMPlexVecGetClosure(dm,section,dsfield->vec,c,&closureSize,&elem);CHKERRQ(ierr); 148 if (B) { 149 if (type == PETSC_SCALAR) { 150 PetscScalar *cB = &((PetscScalar *) B)[nc * nq * i]; 151 152 DMFieldDSdot(cB,fB,elem,nq,feDim,nc,(PetscScalar)); 153 } else { 154 PetscReal *cB = &((PetscReal *) B)[nc * nq * i]; 155 156 DMFieldDSdot(cB,fB,elem,nq,feDim,nc,PetscRealPart); 157 } 158 } 159 if (D) { 160 if (type == PETSC_SCALAR) { 161 PetscScalar *cD = &((PetscScalar *) D)[nc * nq * dim * i]; 162 163 DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),(PetscScalar)); 164 } else { 165 PetscReal *cD = &((PetscReal *) D)[nc * nq * dim * i]; 166 167 DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),PetscRealPart); 168 } 169 } 170 if (H) { 171 if (type == PETSC_SCALAR) { 172 PetscScalar *cH = &((PetscScalar *) H)[nc * nq * dim * dim * i]; 173 174 DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),(PetscScalar)); 175 } else { 176 PetscReal *cH = &((PetscReal *) H)[nc * nq * dim * dim * i]; 177 178 DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),PetscRealPart); 179 } 180 } 181 ierr = DMPlexVecRestoreClosure(dm,section,dsfield->vec,c,&closureSize,&elem);CHKERRQ(ierr); 182 } 183 ierr = PetscFERestoreTabulation(fe,nq,qpoints,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);CHKERRQ(ierr); 184 } else {SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented");} 185 if (!isStride) { 186 ierr = ISRestoreIndices(pointIS,&points);CHKERRQ(ierr); 187 } 188 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 189 PetscFunctionReturn(0); 190 } 191 192 static PetscErrorCode DMFieldGetFEInvariance_DS(DMField field, IS pointIS, PetscBool *isConstant, PetscBool *isAffine, PetscBool *isQuadratic) 193 { 194 DMField_DS *dsfield; 195 PetscObject disc; 196 PetscInt h, imin; 197 PetscClassId id; 198 PetscErrorCode ierr; 199 200 PetscFunctionBegin; 201 dsfield = (DMField_DS *) field->data; 202 ierr = ISGetMinMax(pointIS,&imin,NULL);CHKERRQ(ierr); 203 for (h = 0; h < dsfield->height; h++) { 204 PetscInt hEnd; 205 206 ierr = DMPlexGetHeightStratum(field->dm,h,NULL,&hEnd);CHKERRQ(ierr); 207 if (imin < hEnd) break; 208 } 209 ierr = DMFieldDSGetHeightDisc(field,h,&disc);CHKERRQ(ierr); 210 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 211 if (id == PETSCFE_CLASSID) { 212 PetscFE fe = (PetscFE) disc; 213 PetscInt order, maxOrder; 214 PetscBool tensor = PETSC_FALSE; 215 PetscSpace sp; 216 217 ierr = PetscFEGetBasisSpace(fe, &sp);CHKERRQ(ierr); 218 ierr = PetscSpaceGetOrder(sp,&order);CHKERRQ(ierr); 219 ierr = PetscSpacePolynomialGetTensor(sp,&tensor);CHKERRQ(ierr); 220 if (tensor) { 221 PetscInt dim; 222 223 ierr = DMGetDimension(field->dm,&dim);CHKERRQ(ierr); 224 maxOrder = order * dim; 225 } else { 226 maxOrder = order; 227 } 228 if (isConstant) *isConstant = (maxOrder < 1) ? PETSC_TRUE : PETSC_FALSE; 229 if (isAffine) *isAffine = (maxOrder < 2) ? PETSC_TRUE : PETSC_FALSE; 230 if (isQuadratic) *isQuadratic = (maxOrder < 3) ? PETSC_TRUE : PETSC_FALSE; 231 } 232 PetscFunctionReturn(0); 233 } 234 235 static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad) 236 { 237 PetscInt h, dim, imax, imin; 238 DM dm; 239 DMField_DS *dsfield; 240 PetscObject disc; 241 PetscFE fe; 242 PetscClassId id; 243 PetscErrorCode ierr; 244 245 246 PetscFunctionBegin; 247 dm = field->dm; 248 dsfield = (DMField_DS *) field->data; 249 ierr = ISGetMinMax(pointIS,&imax,&imin);CHKERRQ(ierr); 250 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 251 for (h = 0; h <= dim; h++) { 252 PetscInt hStart, hEnd; 253 254 ierr = DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);CHKERRQ(ierr); 255 if (imin >= hStart && imax < hEnd) break; 256 } 257 *quad = NULL; 258 if (h < dsfield->height) { 259 ierr = DMFieldDSGetHeightDisc(field,h,&disc);CHKERRQ(ierr); 260 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 261 if (id != PETSCFE_CLASSID) PetscFunctionReturn(0); 262 fe = (PetscFE) disc; 263 ierr = PetscFEGetQuadrature(fe,quad);CHKERRQ(ierr); 264 ierr = PetscObjectReference((PetscObject)*quad);CHKERRQ(ierr); 265 } 266 PetscFunctionReturn(0); 267 } 268 269 static PetscErrorCode DMFieldInitialize_DS(DMField field) 270 { 271 PetscFunctionBegin; 272 field->ops->destroy = DMFieldDestroy_DS; 273 field->ops->evaluate = DMFieldEvaluate_DS; 274 field->ops->evaluateFE = DMFieldEvaluateFE_DS; 275 field->ops->getFEInvariance = DMFieldGetFEInvariance_DS; 276 field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS; 277 field->ops->view = DMFieldView_DS; 278 PetscFunctionReturn(0); 279 } 280 281 PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field) 282 { 283 DMField_DS *dsfield; 284 PetscErrorCode ierr; 285 286 PetscFunctionBegin; 287 ierr = PetscNewLog(field,&dsfield);CHKERRQ(ierr); 288 field->data = dsfield; 289 ierr = DMFieldInitialize_DS(field);CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec,DMField *field) 294 { 295 DMField b; 296 DMField_DS *dsfield; 297 PetscObject disc = NULL; 298 PetscBool isContainer = PETSC_FALSE; 299 PetscClassId id = -1; 300 PetscInt numComponents = -1, dsNumFields; 301 PetscSection section; 302 PetscErrorCode ierr; 303 304 PetscFunctionBegin; 305 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 306 ierr = PetscSectionGetFieldComponents(section,fieldNum,&numComponents);CHKERRQ(ierr); 307 ierr = DMGetNumFields(dm,&dsNumFields);CHKERRQ(ierr); 308 if (dsNumFields) {ierr = DMGetField(dm,fieldNum,&disc);CHKERRQ(ierr);} 309 if (disc) { 310 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 311 isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE; 312 } 313 if (!disc || isContainer) { 314 MPI_Comm comm = PetscObjectComm((PetscObject) dm); 315 PetscInt cStart, cEnd, dim; 316 PetscInt localConeSize = 0, coneSize; 317 PetscFE fe; 318 PetscDualSpace Q; 319 PetscSpace P; 320 DM K; 321 PetscQuadrature quad, fquad; 322 PetscBool isSimplex; 323 324 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 325 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 326 if (cEnd > cStart) { 327 ierr = DMPlexGetConeSize(dm, cStart, &localConeSize);CHKERRQ(ierr); 328 } 329 ierr = MPI_Allreduce(&localConeSize,&coneSize,1,MPIU_INT,MPI_MAX,comm);CHKERRQ(ierr); 330 isSimplex = (coneSize == (dim + 1)) ? PETSC_TRUE : PETSC_FALSE; 331 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 332 ierr = PetscSpaceSetOrder(P, 1);CHKERRQ(ierr); 333 ierr = PetscSpaceSetNumComponents(P, numComponents);CHKERRQ(ierr); 334 ierr = PetscSpaceSetType(P,PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 335 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 336 ierr = PetscSpacePolynomialSetTensor(P, isSimplex ? PETSC_FALSE : PETSC_TRUE);CHKERRQ(ierr); 337 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 338 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 339 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 340 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 341 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 342 ierr = DMDestroy(&K);CHKERRQ(ierr); 343 ierr = PetscDualSpaceSetNumComponents(Q, numComponents);CHKERRQ(ierr); 344 ierr = PetscDualSpaceSetOrder(Q, 1);CHKERRQ(ierr); 345 ierr = PetscDualSpaceLagrangeSetTensor(Q, isSimplex ? PETSC_FALSE : PETSC_TRUE);CHKERRQ(ierr); 346 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 347 ierr = PetscFECreate(comm, &fe);CHKERRQ(ierr); 348 ierr = PetscFESetType(fe,PETSCFEBASIC);CHKERRQ(ierr); 349 ierr = PetscFESetBasisSpace(fe, P);CHKERRQ(ierr); 350 ierr = PetscFESetDualSpace(fe, Q);CHKERRQ(ierr); 351 ierr = PetscFESetNumComponents(fe, numComponents);CHKERRQ(ierr); 352 ierr = PetscFESetUp(fe);CHKERRQ(ierr); 353 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 354 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 355 if (isSimplex) { 356 ierr = PetscDTGaussJacobiQuadrature(dim, 1, 1, -1.0, 1.0, &quad);CHKERRQ(ierr); 357 ierr = PetscDTGaussJacobiQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);CHKERRQ(ierr); 358 } 359 else { 360 ierr = PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, &quad);CHKERRQ(ierr); 361 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);CHKERRQ(ierr); 362 } 363 ierr = PetscFESetQuadrature(fe, quad);CHKERRQ(ierr); 364 ierr = PetscFESetFaceQuadrature(fe, fquad);CHKERRQ(ierr); 365 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 366 ierr = PetscQuadratureDestroy(&fquad);CHKERRQ(ierr); 367 disc = (PetscObject) fe; 368 } else { 369 ierr = PetscObjectReference(disc);CHKERRQ(ierr); 370 } 371 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 372 if (id == PETSCFE_CLASSID) { 373 PetscFE fe = (PetscFE) disc; 374 375 ierr = PetscFEGetNumComponents(fe,&numComponents);CHKERRQ(ierr); 376 } else {SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented");} 377 ierr = DMFieldCreate(dm,numComponents,DMFIELD_VERTEX,&b);CHKERRQ(ierr); 378 ierr = DMFieldSetType(b,DMFIELDDS);CHKERRQ(ierr); 379 dsfield = (DMField_DS *) b->data; 380 dsfield->fieldNum = fieldNum; 381 ierr = DMGetDimension(dm,&dsfield->height);CHKERRQ(ierr); 382 dsfield->height++; 383 ierr = PetscCalloc1(dsfield->height,&dsfield->disc);CHKERRQ(ierr); 384 dsfield->disc[0] = disc; 385 ierr = PetscObjectReference((PetscObject)vec);CHKERRQ(ierr); 386 dsfield->vec = vec; 387 *field = b; 388 389 PetscFunctionReturn(0); 390 } 391